# Adaptive Evolution of Virulence-Related Traits in a Susceptible-Infected Model with Treatment

Authors: Xinzhu Meng; Yang Yang; Shengnan Zhao
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

Evolution problem is now a hot topic in the mathematical biology field. This paper investigates the adaptive evolution of pathogen virulence in a susceptible-infected (SI) model under drug treatment. We explore the evolution of a continuous trait, virulence of a pathogen, and consider virulence-dependent cure rate (recovery rate) that dramatically affects the outcome of evolution. With the methods of critical function analysis and adaptive dynamics, we identify the evolutionary conditions for continuously stable strategies, evolutionary repellers, and evolutionary branching points. First, the results show that a high-intensity strength drug treatment can result in evolutionary branching and the evolution of pathogen strains will tend towards a higher virulence with the increase of the strength of the treatment. Second, we use the critical function analysis to investigate the evolution of virulence-related traits and show that evolutionary outcomes strongly depend on the shape of the trade-off between virulence and transmission. Third, after evolutionary branching, the two infective species will evolve to an evolutionarily stable dimorphism at which they can continue to coexist, and no further branching is possible, which is independent of the cure rate function.


---

## Body

## 1. Introduction

Virus leads to the spread of the disease, which attracts the attention and research of many researchers [1–10]. In this paper, we will use a classical susceptible-infected (S-I) epidemic model and assume that all hosts are equally susceptible [1, 2]. Once infected, the host may either die or recover from drug curing, and then, the host becomes fully susceptible again. The model considers a pathogen which increases the mortality of infected hosts by its virulence and whose rate of transmission to uninfected hosts depends on its virulence. High virulence of the pathogen reduces cure rates of infected individuals. There is much good work on the evolution between transmission and virulence [1–5], but the work on the evolutionary adaptive dynamics between drug treatment and virulence is very little (note that this situation is very common in evolutionary ecology and supported by empirical data [6]). Our main aim is to explore how the strength of treatment affects the outcome of evolution. The outcome of evolution often depends on traits which influence infectious ability such as virulence. Therefore, as a phenotype trait, virulence of the pathogen will be taken as an object of study. The evolution of pathogens recently has focused especially on influence virulence and on reasons for the coexistence of parasite strains [2, 7, 8]. Here, we develop an epidemiological model to study the effect of drug treatment on the evolution of pathogen virulence. The model is based upon adaptive dynamics [11–19], in which we explore not only the coexistence of parasite strains but also continuously stable strategy and evolutionary branching.The rest of the paper is organized as follows. In Section2, a susceptible-infected model under drug treatment is first proposed and then the fitness of mutant strains will be deduced. In Section 3, the critical function analysis and adaptive dynamics are carried out to investigate evolutionary dynamics in a monomorphic state. Section 4 discusses the dimorphic coexistence and coevolution of two parasite strains. We then conclude our paper in Section 5 with discussions.


## 2. Model and Demographic Properties

The epidemic model itself for influenza and venereal disease is of the simplest kinds, it is assumed that the cured infective hosts recover and become susceptible, and thus infection is not permanent as they have recovery rate. The infective individuals, although possibly subject to a higher mortality and lower fertility than the susceptible, are otherwise active in the population. The model considers a pathogen which increases the mortality of infected hosts by its virulencex, and an increase in transmission rate can only evolve with a parallel increase in virulence. In this model, we take into account treating, and the recovery becomes susceptible after the infective individuals are cured. It is assumed that an increase in virulence can evolve with a decrease in cure rate (recovery rate).


### 2.1. Monomorphic Evolutionary Dynamics

In this subsection, in case of a monomorphic infective population with traitx, the host dynamics becomes
(1)S′(t)=λN-μ(N)S-β(x)SI+h(x)I,I′(t)=β(x)SI-μ(N)I-α(x)I-h(x)I,
where S and I, respectively, represent the density of the susceptible population and the infected population at time t, with total population density N=S+I,λ is per capita birth rate, and h(x) is the cure rate (the recovery rate) related to virulence of the pathogen. The population is subjected to density-dependent mortality given by the continuous, monotone increasing function μ(N), the infectives suffer also an additional mortality α(x) due to disease, and the transmission rate β(x) depends on virulence of the pathogen.Concerning the effect of virulencex of the pathogen on the transmission rate β(x), mortality α(x) due to disease and the cure rate h(x), biological considerations dictate the properties: β(x) and α(x) are continuous and increasing functions, while h(x) is a continuous and decreasing function because strong virulence means large transmission rate, large mortality due to disease, and small cure rate.In this model we assume linear density-dependent mortality. Letμ(N)=A+BN such that, in absence of the pathogen, the host exhibits logistic growth. Without loss of generality we adopt the values A=B=1. Now setting the right-hand sides of (1) to 0, when
(2)β(x)>μ(N*)+α(x)+h(x)N*
is satisfied, the system (1) has a unique stable positive equilibrium (N*(x),S*(x),I*(x)) [14], where(3)N*(x)=(α(x)-β(x)+λβ(x)-α(x)β(x)+-4(λ+h(x))(1+h(x)+α(x))β(x)2+((1+λ+h(x))β(x)+α(x)(1+β(x)))2h+-4(λ+h(x))(1+h(x)+α(x))β(x)2+((1+λ+h(x))β(x)+α(x)(1+β(x)))2)(2β(x))-1,S*(x)=(α(x)+β(x)+λβ(x)+2h(x)β(x)+α(x)β(x)+-4(λ+h(x))(1+h(x)+α(x))β(x)2+((1+λ+2h(x))β(x)+α(x)(1+β(x)))2h+-4(λ+h(x))(1+h(x)+α(x))β(x)2+((1+λ+2h(x))β(x)+α(x)(1+β(x)))2)(2β(x)2)-1,I*(x)=N*(x)-S*(x).To find the invasion fitness for a mutant infected population, we extend the resident population model (1) by introducing a mutant infected population Iy. If the mutant infected population with a slightly different trait y enters into the resident susceptible-infected system at a low density, the resident-mutant population dynamics is given by
(4)S′(t)=λN-μ(N)S-β(x)SI-β(y)SIy+h(x)I+h(y)Iy,I′(t)=β(x)SI-μ(N)I-α(x)I-h(x)I,Iy′(t)=β(y)SIy-μ(N)Iy-α(y)Iy-h(y)Iy.
Here, Iy denotes the population density of the mutant infected population at time t. One ecological equilibrium of the resident-mutant model (4) is (S*(x),I*(x),0), where S*(x) and I*(x) are the same as those in (3). Just before the occurrence of the mutation, the resident populations are close to the steady state (3) of the resident population model (1) or land at it. Thus, just after the small and rare mutation, the resident and mutant populations are close to the ecological equilibrium of model (4), I*(x) and 0, respectively. However, the stability of the equilibrium (S*(x),I*(x),0) of model (4) determines whether or not the mutant predators can invade. The stability of the extinction equilibrium (S*(x),I*(x),0) is governed by the following Jacobian matrix J, that results from straight forward linearization:(5)J=(λ-∂μ(N*)∂S*-μ(N*)-β(x)I*λ-∂μ(N*)∂I*S*-β(x)S*+h(x)λ-∂μ(N*)∂Iy*S*-β(y)S*+h(y)β(x)I*-∂μ(N*)∂S*I*-∂μ(N*)∂I*I*-∂μ(N*)∂Iy*I*00β(y)S*-μ(N*)-α(y)-h(y)).To show its structure, we rewrite matrixJ in a way that uses its block-triangular form:
(6)J=(JresJ10Jmut),
where(7)Jres=[λ-∂μ(N*)∂S*-μ(N*)-β(x)I*λ-∂μ(N*)∂I*S*-β(x)S*+h(x)β(x)I*-∂μ(N*)∂S*I*-∂μ(N*)∂I*I*],0=(0,0),J1=[λ-∂μ(N*)∂Iy*S*-β(y)S*+h(y),-∂μ(N*)∂Iy*I*]T,Jmut=β(y)S*-μ(N*)-α(y)-h(y).Here J is an upper triangular matrix, so its eigenvalues are given by the eigenvalues of the diagonal blocks, Jres and Jmut. We have already acquired that if (2) holds, (S*(x),I*(x)) is globally asymptotically stable before mutant predators appear, which implies the eigenvalues of Jres have negative real parts. Therefore, if Jmut<0, then all eigenvalues of J have negative real parts, the equilibrium (S*(x),I*(x),0) of system (4) is stable, and the mutant cannot invade. Otherwise, they can invade. So the relevant is the single element of Jmut, which we define as
(8)sx(y)=β(y)S*(x)-μ(N*(x))-α(y)-h(y),
where sx(x)=0.Note thatsx(y) is just the per capita growth rate of hosts infected with the mutant strain. If sx(y)>0, the abundance of mutant infected strains will initially increase, the mutant strains can successfully invade the resident strains; while if sx(y)<0, mutant infected strains cannot invade, and they are doomed to extinct. Thus, sx(y) is, namely, invasion fitness [15], so we think the invasion fitness determines the fate of mutant infected strains.From (8), we can clearly see that there are two aspects of the environment, into which a mutant enters, that influence its invasion fitness: the density of the susceptible and the total density of hosts (which acts through the density-dependent mortality). These are both functions of the resident strain and will thus act as feedback variables in the evolutionary process. When mutations are very small, trait y and trait x are similar to each other and we have as a linear approximation of the mutant’s fitness [13]
(9)sx(y)=sx(x)+D(x)(y-x),
where D(x)=(∂sx(y)/∂y)|y=x.IfD(x)>0, it is clear that only mutants with y>x can invade and take over resident infected strains, whereas if D(x)<0, then it is only possible for mutants with y<x. The virulence thus evolves in the direction of the local fitness gradient until it reaches the neighborhood of a strategy for which D(x) is zero. A strategy x* for which the local fitness gradient is zero is called an evolutionarily singular strategy [13]; that is,
(10)D(x*)=∂sx(y)∂y|y=x=x*=β′(x*)S*(x*)-α′(x*)-h′(x*)=0⟹β′(x*)=α′(x*)+h′(x*)S*(x*).Since mutations are random and very small, the evolutionary model of virulencex can be given by
(11)dxdt=12δσ2I*(x)D(x),
where δ is the probability of individual mutation which is born in resident infected population; ρ2 is the variance of mutation distribution of infected population; (1/2)δρ2 is mutational rate of the infected; I*(x) is the population density of the resident infected at positive equilibrium.A singular strategyx* is evolutionarily stable strategy (ESS) [16] if no nearby mutant can invade, in other words, if sx*(y)<0 for all y≠x* in a neighborhood of x*. If the singular strategy x* is a local fitness maximum, it is a local ESS satisfying
(12)∂2sx(y)∂y2|y=x=x*=β′′(x*)α′(x*)+h′(x*)β′(x*)-α′′(x*)-h′′(x*)<0.
Since at x* the sign of the local fitness gradient changes from positive to negative, D(x) is a decreasing function with respect to x. Directional evolution of resident infected strain can approach singular strategy, so the singularity must be convergence stable [15], if
(13)dD(x)dx|y=x=x*=ddx[∂sx(y)∂y|y=x]|x=x*=∂2sx(y)∂y2|y=x=x*+∂2sx(y)∂y∂x|y=x=x*<0.
A singularity that is not only evolutionarily stable but also convergence stable is said to be a continuously stable strategy (CSS) and will represent a stable endpoint of evolution. In this case, the continuously stable strategy satisfies (12) and (13), and thus the resident pathogen has a maximum fitness. However, if (13) holds but (12) does not, gradual evolution will proceed towards a singularity at which any nearby mutant can invade such that the resident pathogen can undergo evolutionary branching and in which case the singularity will become an evolutionary branching point. The evolving pathogen will turn dimorphic after evolutionary branching. At least initially, the dimorphic pathogens will evolve away from the singularity and thus apart from one another. If (13) does not hold, the singularity is a repeller.In this model, we choose the cure rate function, transmission rate function, and mortality function of infected hosts like this:(14)h(x)=c(1-11+vexp⁡(-x/(1+p))),β(x)=axb+x(1-Kexp⁡(-(x-μ)2σ2)),α(x)=x,
where p represents the strength of treatment disturbance [7, 17].We use “pairwise invasibility plot” (PIP) to study how the strength of treatment affects the outcome of evolution. So we choose three different values ofp, we have plotted three PIPs in Figure 1, in which x-axis and y-axis represent the resident strategy and the mutant strategy, respectively. We look along a vertical line through the evolutionary singular strategy on the x-axis. By the conditions of ESS-stability, we know if the singular strategy x* is ESS-stable, it will satisfy that the vertical line through x* lies completely inside a region marked “−” corresponding to mutants for which sx(y)<0 and which therefore cannot invade; we can see when p=0.5, the singular strategy is an ESS (see Figure 1(a)). Figure 1(a) shows that when p=0.5, the singular strategy is not only ESS-stable but also convergence stable, so the singular strategy x* is a CSS. By the conditions of evolutionary branching, we know if the singular strategy is convergence stable, but not evolutionarily stable, the parts of vertical line inside a region marked “+” correspond to strategies on the y-axis for which sx(y)>0 and hence denote potentially invading mutants. Obviously when p=6 and p=8, the mutant pathogen strains can invade, so evolutionary branching occurs (see Figures 1(b) and 1(c)).Pairwise invasibility plots. The resident and mutant strategies are denoted byx and y, respectively. The shaded areas marked with “+” indicate the combinations of x and y for which the fitness of the mutant, sx(y), is positive. Arrows indicate the direction of disruptive selection which implies convergent stability. Black bullet A = evolutionarily stable singular strategy, red bullet B = evolutionary branching point, and blue bullet C = repeller. (a) Only one singular strategy is a CSS when p=0.5; (b) only one singular strategy is an evolutionary branching point when p=6; and (c) multisingular strategies p=8. Model parameters: a=200, b=10, μ=6, σ=10, c=2, v=0.8, K=0.3, and λ=10.25.
(a)
(b)
(c)


### 2.2. The Effect of Treatment Disturbance Strength on Evolutionary Dynamics

Correspondingly, we plot the bifurcation diagram in Figure2, and through the diagram we can see the singular trait values (the virulence x) as a function of the strength of the treatment p, which shows the stability properties of the singular strategy will change due to the changing values of p. An increase in values of p will parallel increase in the singular trait values (the virulence x). This shows an increase in the strength of the treatment can lead to stable evolution towards larger virulence x. This means the pathogen strain will have immunity against drugs. Then the strain outcompetes its rivals and gives birth to the higher virulence mutant pathogen strains. In the lower part, when p is near zero, that is to say, the cure rate is very small, then there is an ESS, but it loses ESS-stability when p>1.49; then all convergence stable singularities are branching points. It is easy to see that a low strength drug treatment may induce a continuously stable strategy, while a high strength treatment can induce evolutionary branching. Moreover, there is a complicated evolutionary dynamics including evolutionary branching, CSS, and repeller when p>7.5.Figure 2
Bifurcation diagram for evolutionarily singular strategyx and recovery disturbance strength p. Red solid lines indicate unstable singular strategy x (i.e., evolutionary branching point), black solid lines indicate the CSS, and red dashed line indicates the repeller. The parameter values are a=200, b=10, μ=6, σ=10, c=2, v=0.8, K=0.3, and λ=10.25.


## 2.1. Monomorphic Evolutionary Dynamics

In this subsection, in case of a monomorphic infective population with traitx, the host dynamics becomes
(1)S′(t)=λN-μ(N)S-β(x)SI+h(x)I,I′(t)=β(x)SI-μ(N)I-α(x)I-h(x)I,
where S and I, respectively, represent the density of the susceptible population and the infected population at time t, with total population density N=S+I,λ is per capita birth rate, and h(x) is the cure rate (the recovery rate) related to virulence of the pathogen. The population is subjected to density-dependent mortality given by the continuous, monotone increasing function μ(N), the infectives suffer also an additional mortality α(x) due to disease, and the transmission rate β(x) depends on virulence of the pathogen.Concerning the effect of virulencex of the pathogen on the transmission rate β(x), mortality α(x) due to disease and the cure rate h(x), biological considerations dictate the properties: β(x) and α(x) are continuous and increasing functions, while h(x) is a continuous and decreasing function because strong virulence means large transmission rate, large mortality due to disease, and small cure rate.In this model we assume linear density-dependent mortality. Letμ(N)=A+BN such that, in absence of the pathogen, the host exhibits logistic growth. Without loss of generality we adopt the values A=B=1. Now setting the right-hand sides of (1) to 0, when
(2)β(x)>μ(N*)+α(x)+h(x)N*
is satisfied, the system (1) has a unique stable positive equilibrium (N*(x),S*(x),I*(x)) [14], where(3)N*(x)=(α(x)-β(x)+λβ(x)-α(x)β(x)+-4(λ+h(x))(1+h(x)+α(x))β(x)2+((1+λ+h(x))β(x)+α(x)(1+β(x)))2h+-4(λ+h(x))(1+h(x)+α(x))β(x)2+((1+λ+h(x))β(x)+α(x)(1+β(x)))2)(2β(x))-1,S*(x)=(α(x)+β(x)+λβ(x)+2h(x)β(x)+α(x)β(x)+-4(λ+h(x))(1+h(x)+α(x))β(x)2+((1+λ+2h(x))β(x)+α(x)(1+β(x)))2h+-4(λ+h(x))(1+h(x)+α(x))β(x)2+((1+λ+2h(x))β(x)+α(x)(1+β(x)))2)(2β(x)2)-1,I*(x)=N*(x)-S*(x).To find the invasion fitness for a mutant infected population, we extend the resident population model (1) by introducing a mutant infected population Iy. If the mutant infected population with a slightly different trait y enters into the resident susceptible-infected system at a low density, the resident-mutant population dynamics is given by
(4)S′(t)=λN-μ(N)S-β(x)SI-β(y)SIy+h(x)I+h(y)Iy,I′(t)=β(x)SI-μ(N)I-α(x)I-h(x)I,Iy′(t)=β(y)SIy-μ(N)Iy-α(y)Iy-h(y)Iy.
Here, Iy denotes the population density of the mutant infected population at time t. One ecological equilibrium of the resident-mutant model (4) is (S*(x),I*(x),0), where S*(x) and I*(x) are the same as those in (3). Just before the occurrence of the mutation, the resident populations are close to the steady state (3) of the resident population model (1) or land at it. Thus, just after the small and rare mutation, the resident and mutant populations are close to the ecological equilibrium of model (4), I*(x) and 0, respectively. However, the stability of the equilibrium (S*(x),I*(x),0) of model (4) determines whether or not the mutant predators can invade. The stability of the extinction equilibrium (S*(x),I*(x),0) is governed by the following Jacobian matrix J, that results from straight forward linearization:(5)J=(λ-∂μ(N*)∂S*-μ(N*)-β(x)I*λ-∂μ(N*)∂I*S*-β(x)S*+h(x)λ-∂μ(N*)∂Iy*S*-β(y)S*+h(y)β(x)I*-∂μ(N*)∂S*I*-∂μ(N*)∂I*I*-∂μ(N*)∂Iy*I*00β(y)S*-μ(N*)-α(y)-h(y)).To show its structure, we rewrite matrixJ in a way that uses its block-triangular form:
(6)J=(JresJ10Jmut),
where(7)Jres=[λ-∂μ(N*)∂S*-μ(N*)-β(x)I*λ-∂μ(N*)∂I*S*-β(x)S*+h(x)β(x)I*-∂μ(N*)∂S*I*-∂μ(N*)∂I*I*],0=(0,0),J1=[λ-∂μ(N*)∂Iy*S*-β(y)S*+h(y),-∂μ(N*)∂Iy*I*]T,Jmut=β(y)S*-μ(N*)-α(y)-h(y).Here J is an upper triangular matrix, so its eigenvalues are given by the eigenvalues of the diagonal blocks, Jres and Jmut. We have already acquired that if (2) holds, (S*(x),I*(x)) is globally asymptotically stable before mutant predators appear, which implies the eigenvalues of Jres have negative real parts. Therefore, if Jmut<0, then all eigenvalues of J have negative real parts, the equilibrium (S*(x),I*(x),0) of system (4) is stable, and the mutant cannot invade. Otherwise, they can invade. So the relevant is the single element of Jmut, which we define as
(8)sx(y)=β(y)S*(x)-μ(N*(x))-α(y)-h(y),
where sx(x)=0.Note thatsx(y) is just the per capita growth rate of hosts infected with the mutant strain. If sx(y)>0, the abundance of mutant infected strains will initially increase, the mutant strains can successfully invade the resident strains; while if sx(y)<0, mutant infected strains cannot invade, and they are doomed to extinct. Thus, sx(y) is, namely, invasion fitness [15], so we think the invasion fitness determines the fate of mutant infected strains.From (8), we can clearly see that there are two aspects of the environment, into which a mutant enters, that influence its invasion fitness: the density of the susceptible and the total density of hosts (which acts through the density-dependent mortality). These are both functions of the resident strain and will thus act as feedback variables in the evolutionary process. When mutations are very small, trait y and trait x are similar to each other and we have as a linear approximation of the mutant’s fitness [13]
(9)sx(y)=sx(x)+D(x)(y-x),
where D(x)=(∂sx(y)/∂y)|y=x.IfD(x)>0, it is clear that only mutants with y>x can invade and take over resident infected strains, whereas if D(x)<0, then it is only possible for mutants with y<x. The virulence thus evolves in the direction of the local fitness gradient until it reaches the neighborhood of a strategy for which D(x) is zero. A strategy x* for which the local fitness gradient is zero is called an evolutionarily singular strategy [13]; that is,
(10)D(x*)=∂sx(y)∂y|y=x=x*=β′(x*)S*(x*)-α′(x*)-h′(x*)=0⟹β′(x*)=α′(x*)+h′(x*)S*(x*).Since mutations are random and very small, the evolutionary model of virulencex can be given by
(11)dxdt=12δσ2I*(x)D(x),
where δ is the probability of individual mutation which is born in resident infected population; ρ2 is the variance of mutation distribution of infected population; (1/2)δρ2 is mutational rate of the infected; I*(x) is the population density of the resident infected at positive equilibrium.A singular strategyx* is evolutionarily stable strategy (ESS) [16] if no nearby mutant can invade, in other words, if sx*(y)<0 for all y≠x* in a neighborhood of x*. If the singular strategy x* is a local fitness maximum, it is a local ESS satisfying
(12)∂2sx(y)∂y2|y=x=x*=β′′(x*)α′(x*)+h′(x*)β′(x*)-α′′(x*)-h′′(x*)<0.
Since at x* the sign of the local fitness gradient changes from positive to negative, D(x) is a decreasing function with respect to x. Directional evolution of resident infected strain can approach singular strategy, so the singularity must be convergence stable [15], if
(13)dD(x)dx|y=x=x*=ddx[∂sx(y)∂y|y=x]|x=x*=∂2sx(y)∂y2|y=x=x*+∂2sx(y)∂y∂x|y=x=x*<0.
A singularity that is not only evolutionarily stable but also convergence stable is said to be a continuously stable strategy (CSS) and will represent a stable endpoint of evolution. In this case, the continuously stable strategy satisfies (12) and (13), and thus the resident pathogen has a maximum fitness. However, if (13) holds but (12) does not, gradual evolution will proceed towards a singularity at which any nearby mutant can invade such that the resident pathogen can undergo evolutionary branching and in which case the singularity will become an evolutionary branching point. The evolving pathogen will turn dimorphic after evolutionary branching. At least initially, the dimorphic pathogens will evolve away from the singularity and thus apart from one another. If (13) does not hold, the singularity is a repeller.In this model, we choose the cure rate function, transmission rate function, and mortality function of infected hosts like this:(14)h(x)=c(1-11+vexp⁡(-x/(1+p))),β(x)=axb+x(1-Kexp⁡(-(x-μ)2σ2)),α(x)=x,
where p represents the strength of treatment disturbance [7, 17].We use “pairwise invasibility plot” (PIP) to study how the strength of treatment affects the outcome of evolution. So we choose three different values ofp, we have plotted three PIPs in Figure 1, in which x-axis and y-axis represent the resident strategy and the mutant strategy, respectively. We look along a vertical line through the evolutionary singular strategy on the x-axis. By the conditions of ESS-stability, we know if the singular strategy x* is ESS-stable, it will satisfy that the vertical line through x* lies completely inside a region marked “−” corresponding to mutants for which sx(y)<0 and which therefore cannot invade; we can see when p=0.5, the singular strategy is an ESS (see Figure 1(a)). Figure 1(a) shows that when p=0.5, the singular strategy is not only ESS-stable but also convergence stable, so the singular strategy x* is a CSS. By the conditions of evolutionary branching, we know if the singular strategy is convergence stable, but not evolutionarily stable, the parts of vertical line inside a region marked “+” correspond to strategies on the y-axis for which sx(y)>0 and hence denote potentially invading mutants. Obviously when p=6 and p=8, the mutant pathogen strains can invade, so evolutionary branching occurs (see Figures 1(b) and 1(c)).Pairwise invasibility plots. The resident and mutant strategies are denoted byx and y, respectively. The shaded areas marked with “+” indicate the combinations of x and y for which the fitness of the mutant, sx(y), is positive. Arrows indicate the direction of disruptive selection which implies convergent stability. Black bullet A = evolutionarily stable singular strategy, red bullet B = evolutionary branching point, and blue bullet C = repeller. (a) Only one singular strategy is a CSS when p=0.5; (b) only one singular strategy is an evolutionary branching point when p=6; and (c) multisingular strategies p=8. Model parameters: a=200, b=10, μ=6, σ=10, c=2, v=0.8, K=0.3, and λ=10.25.
(a)
(b)
(c)


## 2.2. The Effect of Treatment Disturbance Strength on Evolutionary Dynamics

Correspondingly, we plot the bifurcation diagram in Figure2, and through the diagram we can see the singular trait values (the virulence x) as a function of the strength of the treatment p, which shows the stability properties of the singular strategy will change due to the changing values of p. An increase in values of p will parallel increase in the singular trait values (the virulence x). This shows an increase in the strength of the treatment can lead to stable evolution towards larger virulence x. This means the pathogen strain will have immunity against drugs. Then the strain outcompetes its rivals and gives birth to the higher virulence mutant pathogen strains. In the lower part, when p is near zero, that is to say, the cure rate is very small, then there is an ESS, but it loses ESS-stability when p>1.49; then all convergence stable singularities are branching points. It is easy to see that a low strength drug treatment may induce a continuously stable strategy, while a high strength treatment can induce evolutionary branching. Moreover, there is a complicated evolutionary dynamics including evolutionary branching, CSS, and repeller when p>7.5.Figure 2
Bifurcation diagram for evolutionarily singular strategyx and recovery disturbance strength p. Red solid lines indicate unstable singular strategy x (i.e., evolutionary branching point), black solid lines indicate the CSS, and red dashed line indicates the repeller. The parameter values are a=200, b=10, μ=6, σ=10, c=2, v=0.8, K=0.3, and λ=10.25.


## 3. Critical Function and Evolutionary Analysis

As stated above, we specify theα(x) and h(x) as (14), the parameters v<1,c and p are constants, and the values of K decide the convexity of the trade-off function β(x); we know if the trade-off function β(x) is concave, the evolutionary dynamics of this model leads to ESS. However, this may not be the case when the trade-off is partly convex, for which such a singularity can constitute an evolutionary branching point or repeller [7]. From the second equation of system (1), we have
(15)S*(x)=μ(N*(x))+α(x)+h(x)β(x).As such (10) constitutes a differential equation which together with (15) can be expressed as
(16)βcrit′(x)=βcrit(x)(α′(x)+h′(x))μ(N*(x))+α(x)+h(x).The solutions of (16) define curves in the x-βcrit-plane that have the slopes required by (10) at all points, with different initial values yielding alternative curves. These curves are referred to as critical functions [7, 8], and singular strategies are those values of virulence for which the trade-off is tangential to any such critical function. Then we superpose a trade-off function β(x) onto the critical functions, which confirms the existence and identifies the directions of the singularity for this and any other choice of trade-off functions. When the singularity identified is located at a point where the evolutionary stability can be directly inferred from the curvature of the trade-off at the point of tangency. Thus we have the following theorems.Theorem 1.
If the trade-off function is concave, we superpose such a trade-off functionβ(x) onto the convex critical function; the evolutionary dynamics of this model leads to ESS.Proof.
From (8), we calculate the second derivative of sx(y) with respect to y,
(17)∂2sx(y)∂y2|y=x=x*=β′′(x*)S*(x*)-α′′(x*)-h′′(x*);h(x) is given from (14), and the parameters v<1, so,
(18)h′′(x)=-cv(1+p)2·exp⁡⁡(-x/(1+p))(vexp⁡(-x/(1+p))-1)(1+vexp⁡(-x/(1+p)))3>0.
When the trade-off function β(x) is concave, β′′(x)<0. So, we have
(19)∂2sx(y)∂y2|y=x=x*<0.
Clearly, the conclusion holds up.Theorem 2.
If the trade-off function is convex and at the same time less convex than the critical function at the point of tangency, then the singularity is an evolutionary branching point.Proof.
The trade-off function is convex, which will not satisfy the ESS condition (12). We write out the convergence stable condition (13),
(20)β′′(x*)s*(x*)-α′′(x*)-h′′(x*)+β′(x*)(s*(x*))′<0;
that is,
(21)β′′(x*)<α′′(x*)+h′′(x*)s*(x*)-β′(x*)(s*(x*))′s*(x*).
From (10), we know
(22)β′(x*)=α′(x*)+h′(x*)s*(x*).
Thus we have
(23)β′′(x*)<α′′(x*)+h′′(x*)s*(x*)-(α′(x*)+h′(x*))(s*(x*))′(s*(x*))2,
where on the right hand side is exactly the derivative of (16), that is, βcrit′′(x*). That is, β′′(x*)<βcrit′′(x*). So it can be seen directly from the critical function that if the trade-off function is convex and at the same time less convex than the critical functions at the point of tangency, the singularity at tangency is convergence stable but not evolutionarily stable; that is, it is an evolutionary branching point.In Figure3(a), when K=0.2, the trade-off function is concave, so the singular strategy is an ESS. We can see in Figure 3(b), when K=0.3, the trade-off function is convex and at the same time less convex than the critical function at the point of tangency x*, such a trade-off results in convergence stable but invadable singularity, that is, evolutionary branching point. Critical functions are always convex, it is always possible to find a trade-off function that is locally convex and less convex than the critical function, satisfying the condition for evolutionary branching. The singular strategy is a branching point. Of course this result is consistent with bifurcation diagram in Figure 2.Critical function analysis (a)-(b), (a)K=0.2, (b) K=0.3; black solid lines indicate the CSS, red solid line indicates convergence stable but not evolutionarily stable singular strategy x* (i.e., evolutionary branching point), and parameter values are a=200,b=10,μ=6,σ=10,c=2,v=0.8,λ=10.25, and p=6.
(a)
(b)


## 4. Dimorphic Coexistence

In this section we investigate whether further evolutionary branching or evolutionarily stable coexistence of the two resident strains occurs. After branching, population dynamics of two residents with phenotypic traitsx1 and x2 and a mutant with trait y appears, where the two resident strategies soon grow too far apart near the branching point. Two classes of infected individuals, I1, I2 in a host population represent a pathogen comprising two different strains, and we now proceed to investigate the evolutionary endpoint of such a dimorphic pathogen. In the case of a dimorphic pathogen, I1 and I2, one with strain x1 and the other with strain x2, the population dynamics is given by
(24)S′(t)=λN-μ(N)S-β(x1)SI1-β(x2)SI2+h(x1)I1+h(x2)I2,I1′(t)=β(x1)SI1-μ(N)I1-α(x1)I1-h(x1)I1,I2′(t)=β(x2)SI2-μ(N)I2-α(x2)I2-h(x2)I2.
Andreasen and Pugliese [1] proved that a necessary and sufficient condition for two pathogenic strains x1 and x2 to coexist at a unique and stable equilibrium of (24) is that both strains can successfully invade one another. This equilibrium is given by
(25)S*=α(x1)-α(x2)+h(x1)-h(x2)β(x1)-β(x2),μ(N*)=(α(x1)β(x2)-α(x2)β(x1)-h(x1)β(x2)h-h(x2)β(x1))(β(x1)-β(x2))-1,I1*=(λ+μ(N*)-α(x2))N*+α(x2)S*α(x1)-α(x2),S*+I1*+I2*=N*.
By the same arguments as the above (8), when mutant pathogen with a slightly different strain y enters into the resident susceptible-infected system at a low density, the invasion fitness is given by
(26)sx1,x2(y)=β(y)S*(x1,x2)-μ(N*(x1,x2))-α(y)-h(y).The selection gradientsg1(x1,x2) and g2(x1,x2) are given by
(27)g1(x1,x2)=∂sx1,x2(y)∂y|y=x1=β′(x1)S*(x1,x2)-α′(x1)-h′(x1),(28)g2(x1,x2)=∂sx1,x2(y)∂y|y=x2=β′(x2)S*(x1,x2)-α′(x2)-h′(x2).
Dimorphic evolutionary singularities are coalitions of resident strategies where the fitness gradient of each strategy is zero. That is, any solution (x1*,x2*) of gi(x1*,x2*)=0 for i=1,2 is called evolutionarily singular coalition.When individual mutations are random and sufficiently small, the evolutionary model of traitsx1 and x2 can be given by
(29)dx1dt=12δ1ρ12I1*(x1,x2)g1(x1,x2),dx2dt=12δ2ρ22I2*(x1,x2)g2(x1,x2),
where Ii*(x1,x2)(i=1,2) is the equilibrium population density of the resident infected Ii(i=1,2); (1/2)δiρi2(i=1,2) is mutational rate of the infected Ii(i=1,2).Theorem 3.
Further evolutionary branching cannot happen in this model whether the cure rate function is concave or convex.Proof.
A singular coalition is evolutionarily stable if all its constituent strategies are ESS, that is, if(30)∂2sx1,x2(y)∂y2|y=xi*,x1=x1*,x2=x2*=β′′(xi*)S*(x1*,x2*)-α′′(xi*)-h′′(xi*)<0
for i=1 and 2.
The condition for absolute convergence stability is given by(31)ddxi[∂sx1,x2(y)∂y|y=xi]|x1=x1*,x2=x2*=∂2sx1,x2(y)∂y2|y=xi*,x1=x1*,x2=x2*+∂sx1,x2(y)∂y∂xi|y=xi*,x1=x1*,x2=x2*<0.
First, we calculate the mixed derivatives of invasion fitness at the singular coalition:
(32)∂2sx1,x2(y)∂y∂xi=β′(y)∂S*(x1,x2)∂xi
for i=1 and 2. From (25), we have
(33)∂S*(x1,x2)∂x1|x1=x1*,x2=x2*=((α′(x1*)+h′(x1*))(β(x1*)-β(x2*))h-(α(x1*)-α(x2*)+h(x1*)-h(x2*))β′(x1*))×((β(x1*)-β(x2*))2)-1,(34)∂S*(x1,x2)∂x2|x1=x1*,x2=x2*=((α′(x2*)+h′(x2*))(β(x2*)-β(x1*))h-(α(x2*)-α(x1*)+h(x2*)-h(x1*))β′(x2*))×((β(x2*)-β(x1*))2)-1.
Since the singularity condition can be obtained from (27), we have
(35)β′(x1*)=α′(x1*)+h′(x1*)S*(x1*,x2*)=(α′(x1*)+h′(x1*))(β(x1*)-β(x2*))α(x1*)-α(x2*)+h(x1*)-h(x2*),(36)β′(x2*)=α′(x2*)+h′(x2*)S*(x1*,x2*)=(α′(x2*)+h′(x2*))(β(x2*)-β(x1*))α(x2*)-α(x1*)+h(x2*)-h(x1*).
Take (35) and (36) into (33) and (34), respectively, we have
(37)∂S*(x1,x2)∂x1|x1=x1*,x2=x2*=0,∂S*(x1,x2)∂x2|x1=x1*,x2=x2*=0.
Then the mixed derivatives of invasion fitness can readily be verified to be zero. Thus (31) is equivalent to (30). This shows that further evolutionary branching cannot happen in this model.Now, the dynamics of evolution in dimorphic pathogen is shown through the coexistence plots and evolutionary branching trees. The dynamics of evolution as predicted by the model are conformed by numerical simulations (Figure4). In these simulations, we use three different values of p (as indicated above in each graph). In the case of p=0.5, from Section 2, the PIP in the Figure 1(a) shows that the singular strategy x* is an ESS. At the singularity neither mutants from below nor above can invade and the singularity is therefore also convergence stable, and this trade-off thus supports a CSS. Therefore, x* is an evolutionary endpoint and is not invadable by any other mutant strains. In this case an initially monomorphic pathogen population can stay monomorphic. In the corresponding coexistence plot one can see that no singular dimorphism exists, and evolution will drive any dimorphism back into a monomorphic state [2, 7]. This outcome is illustrated in Figures 4(a) and 4(d).The coexistence plots are depicted, two strains are presented on separate axes and the gray areas define combinations of strains that are mutually invasible and the corresponding evolutionary branching. Shaded areas indicating protected dimorphism are separated by stable (black) and unstable (red) isoclines at which selection gradient vanishes in eitherx1-direction or x2-direction. Red bullet B = initial branching point, black bullet A = evolutionarily stable singularity, black bullet Ei(x1*,x2*) = stable singular coalition, and i=1,2. “•” P = repeller point. Correspondingly, below the coexistence plots we simulate evolutionary trees with branching patterns. (a) and (d) p=0.5; (b) and (e) p=6; (c) and (f) p=8. Other parameter values are c=2, v=0.8, λ=10.25, k=0.3.
(a)
(b)
(c)
(d)
(e)
(f)With the effect of increasingp the singular point x* continues to be convergence stable but not an ESS, and in the vicinity of x* a dimorphism can appear, and such x* is an evolutionary branching point (Figures 1(b) and 1(c)). In this case an initially monomorphic population will be close to the singular branching point and will undergo disruptive selection, then becoming dimorphic and comprising two closely related resident phenotypic traits (the virulences x1* and x2*; see numerical simulations in Figures 4(b) and 4(c)). Finally, the two resident pathogen groups are stable coexistence and completely host specialized. In this case no mutant with nonzero transmission on both hosts is able to invade. In the coexistence plot,the singularity loses its evolutionary stability; CSS turns into a branching point. With a further increase in p, when the values of p are 8, then the shaded areas intersect with main diagonal at three points corresponding to three different evolutionarily singular points which present complex dynamics such as ESS, repeller, and evolutionary branching point (Figure 1(c)). Theorem 3 shows that the population evolves towards a stable dimorphism. As is shown in Figure 4, from the corresponding evolutionary branching we can also obtain this result.


## 5. Discussion

In this paper, we employ the susceptible-infected epidemic model under drug treatment to explore the adaptive evolution of pathogen virulence. Our results show that the strength of treatment noticeably affects the outcome of evolution: continuously stable strategies, evolutionary repellers, and evolutionary branching points. With the methods of critical function analysis and adaptive dynamics, we have shown that evolutionary branching of virulence is possible in a susceptible-infected model with cure rate. Conditions ensuring evolutionary stability and branching for singular strategies and coalitions are derived in both monomorphic and dimorphic environments.We obtain the biological interpretation for the evolution of pathogen virulence that if there is no treatment or small treatment, then the pathogen virulence can be smaller and evolutionarily stable, while a powerful treatment could accelerate pathogen evolution towards higher virulence.First, by using the methods of “pairwise invasibility plot” and critical function analysis, we show that if the trade-off function is concave, the evolutionary dynamics of this model can lead to evolutionarily stable, while if the trade-off is convex then the singularity can be an evolutionary branching point. Second, increase in the strength of the drug treatment can increase evidently in the singular strategy value and eventually turns continuously stable strategy into the branching point. This shows that, with increasing in cure rates, pathogen populations will tend towards the evolutionary branching points with higher virulence. Note that this situation is supported by empirical data [6]. Third, our analysis reveals that further evolutionary branching cannot happen in this model whether the cure rate function is concave or convex. For higher cure rates, after branching has occurred in virulence phenotype, the mutant pathogen becomes a resident pathogen and the two resident pathogen strains stably coexist.



---

*Source: 891401-2014-08-07.xml*