Coevolutionary dynamics of host-pathogen interaction with density-dependent mortality.

This study explores the coevolutionary dynamics of host-pathogen interaction based on a susceptible-infected population model with density-dependent mortality. We assume that both the host's resistance and the pathogen's virulence will adaptively evolve, but there are inevitable costs in terms of host birth rate and disease-related mortality rate. Particularly, it is assumed that both the host resistance and pathogen virulence can affect the transmission rate. By using the approach of adaptive dynamics and numerical simulation, we find that the finally coevolutionary outcome depends on the strength of host-pathogen asymmetric interaction, the curvature of trade-off functions, and the intensity of density-dependent natural mortality. To be specific, firstly, we find that if the strengths of host-pathogen asymmetric interaction and disease-related mortality are relatively weak, or the density-dependent natural mortality is relatively strong, then the host resistance and pathogen virulence will evolve to a continuously stable strategy. However, if the strength of host-pathogen asymmetric interaction and disease-related mortality becomes stronger, then the host resistance and pathogen virulence will evolve periodically. Secondly, we find that if the intensities of both the birth rate trade-off function and the density-dependent natural mortality are relatively weak, but the strength of host-pathogen asymmetric interaction becomes relatively strong, then the evolution of host resistance will have a relatively strongly accelerating benefit, the evolutionary branching of host resistance will first arise. However, if the strength of host-pathogen asymmetric interaction is relatively weak, but the intensity of the trade-off function of disease-related mortality becomes relatively strong, then the evolution of pathogen virulence will have a relatively strongly decelerating cost, and the evolutionary branching of pathogen virulence will first arise. Thirdly, after the evolutionary branching of host resistance and pathogen virulence, we further study the coevolutionary dynamics of two-hosts-one-pathogen interaction and one-host-two-pathogens interaction. We find that if the evolutionary branching of host resistance arises firstly, then the finally evolutionary outcome contains a dimorphic host and a monomorphic pathogen population. If the evolutionary branching of pathogen virulence arises firstly, then the finally evolutionary outcome may contain a monomorphic host and a dimorphic pathogen population.

Introduction

The interaction between host and pathogen often has large effect on individual fitness and can significantly alter the evolutionary trajectory of a specie (Tibayrenc 2011). However, we lack an understanding of the evolutionary process that generates and maintains the trait diversity. Besides, the long-term behavior of host-pathogen interaction is likely to depend on the interplay of both species evolutionary traits, that is, the coevolutionary dynamics (Best et al. 2009). Therefore, the coevolutionary dynamics need to be considered more generally, which has important implication for our understanding the coevolutionary mechanism of host-pathogen interaction.

The importance of host-pathogen interaction has led to large theoretical studies that explored the evolution of a single phenotypic trait, pathogen virulence (May and Anderson 1983; Geritz et al. 1998; Tompkins et al. 2002; Day and Proulx 2004; Alizon et al. 2009; Boldin and Diekmann 2008; Duffy and Sivars-Becker 2007; Bowers et al. 2005) or host resistance (Anderson and May 1982; Antonovics and Thrall 1994; Bonds 2006; Boots et al. 2009; Boots and Bowers 1999, 2004; Boots and Haraguchi 1999; Best et al. 2015, 2009; Svennungsen and Kisdi 2009; Gascuel et al. 2013; Kraaijeveld et al. 2012; Mealor and Boots 2006). These studies have established a broad framework for theoretical study and application. Modern theoretical studies on the evolution of pathogen, in which virulence is assumed to be a trade-off against transmission, the so-called ‘virulence-transmission trade-off’, the trade-off theory assumes that a pathogen cannot simultaneously increase transmission and prolong infection, so pathogens are attempting to maximize subject to these costs (May and Anderson 1983; Best et al. 2009). Based on the theory of adaptive dynamics, Alizon et al. (2009) found that in the susceptible-infected (SI) model with density-dependent infection, as long as the trade-off function between transmission rate and pathogen virulence is locally convex, then the evolutionary branching may appear in the locally convex region of the trade-off function. Similar studies have found that in order to obtain the evolutionary branching, the second derivative of the trade-off function near the singular strategy must be less than zero (Boldin and Diekmann 2008; Day and Proulx 2004; Duffy and Sivars-Becker 2007). After branching, the two pathogens with different virulence may coexist. Particularly, when the death rate of host population is density-dependent, but the birth rate is not density-dependent, then the two pathogens with different virulence may coexist (Bowers et al. 2005; Day and Proulx 2004; Geritz et al. 1998; Pugliese 2002), but it is found that the coexistence of two pathogens with different virulence on the population timescale does not guarantee the evolutionary stable coexistence on the long-term evolutionary timescale, one of the branches may evolve to extinction.

Moreover, we have noted that in order to avoid infection by pathogen, the host often needs to change the resistance ability, but this usually requires a cost, the host with higher resistance pays a cost on other life-history traits, such as fertility. This is also supported by theoretical argument and empirical evidence (Anderson and May 1982; Antonovics and Thrall 1994; Bonds 2006; Boots et al. 2009; Boots and Bowers 1999, 2004; Boots and Haraguchi 1999; Best et al. 2015, 2009; Svennungsen and Kisdi 2009; Gascuel et al. 2013; Kraaijeveld et al. 2012; Mealor and Boots 2006). Over the last few decades, much work has been dedicated to understanding the evolutionary mechanism of host resistance, and many different modeling approaches have been developed, such as quantitative genetic method, locus-based method, and adaptive dynamics approach (Best et al. 2015, 2009; Bonds 2006; Boots et al. 2009; Boots and Bowers 1999, 2004; Boots and Haraguchi 1999; Kisdi and Geritz 2016; Landi et al. 2013).

In particular, the approach of adaptive dynamics allows us to examine the evolutionary factors that result in different types of host resistance. Based on the approach of adaptive dynamics (Dieckmann and Law 1996; Geritz et al. 1998; Metz et al. 1992), Antonovics and Thrall (1994) and Bowers et al. (1994) studied the SI model of two-strains, they found a key result of the coexistence of a highly susceptible strain and a highly resistant strain. Boots and Haraguchi (1999) proposed a multi-strain model to study the evolution of host resistance, and they hypothesized that different host strains have different susceptibility to infection strains and that strains with stronger resistance have a lower intrinsic growth rate. The adaptive dynamics approach and pairwise invasibility plot analysis showed that the evolutionary outcome depends on the trade-off function between resistance and the cost. When resistance has a deceleration cost, they found that the evolutionary branch of host resistance may occur.

However, in nature, the interaction between the host and pathogen is usually a coevolutionary process Best et al. (2009), and both the pathogen virulence and host resistance may adaptively evolve. Understanding the coevolutionary mechanism of host-pathogen interaction is one of the key challenges for evolutionary biology. Mathematical models have therefore played an important role in shaping our understanding of host-pathogen coevolution. Recently, many theoretical models are constructed to study the coevolutionary theory of host-pathogen interaction, and find that the population dynamics and the genetic basis of infection have a significant impact on the outcome of host-pathogen coevolution (Buckingham and Ashby 2022). Baalen (1998) and Restif and Koella (2003) investigated the so-called coevolutionary stable state, for which once attained cannot be invaded by other strains. Best et al. (2009) and Ashby et al. (2019) introduced the ecological feedback into the host-pathogen coevolution model to evaluate the impact of ecological evolution itself. Boots et al. (2014) showed that when there is a cost in host resistance and pathogen transmission, epidemiological feedback may produce diversity, but this situation is limited to dimorphism in a wide range of realistic infection scenarios. However, mathematical models that examine the coevolutionary dynamics of quantitative traits in host and pathogen are relatively fewer, especially based on the susceptible-infected population model and the theory of adaptive dynamics (Boots et al. 2009; Best et al. 2010, 2011; Lopez Pascua et al. 2014; Boots et al. 2014; McLeod and Day 2015; Hesse et al. 2015; Best 2018; Alizon 2021). Based on the susceptible-infected population dynamics, Best et al. (2009) showed that highly virulent parasites may evolve due to the coevolutionary process. Best et al. (2010) examined the evolution of host defense to the sterilizing effects of parasites. Best et al. (2011) and Hesse et al. (2015) considered the effect of spatial structure on evolution. Lopez Pascua et al. (2014) explored how the environment influences coevolutionary dynamics. McLeod and Day (2015) considered the implication of avoidance plasticity for host-pathogen coevolution. Alizon (2021) discussed the treatment symptomatic infection and the coevolution of virulence and drug resistance. Best (2018) assessed the impact of predation on the coevolution of costly host resistance and pathogen transmission. Boots et al. (2014) examined how specificity and epidemiology drive the coevolution of static trait diversity in host and parasite. Furthermore, density-dependent mortality has been identified as a key factor in evolutionary dynamics (Andreasen and Pugliese 1995; Zu et al. 2020). Andreasen and Pugliese (1995) discussed that two pathogens can coexist if hosts are subject to density-dependent mortality. Zu et al. (2020) assumed that both susceptible hosts and infected hosts are affected by density-dependent mortality, but only the resistance-related trait of susceptible host can adaptively evolve. They studied the evolutionary mechanism of host resistance to pathogen infection, and found that without density-dependent mortality, the evolutionary branching in the host resistance may not occur. However, in the previous studies, they did not consider the coevolutionary dynamics of host-pathogen interaction when a density-dependent mortality affects all hosts. In addition, after the evolutionary branching of pathogen virulence or host resistance, many previous studies did not continue to explore whether the evolutionary diversity of pathogen virulence or host resistance can be maintained or not for a long time.

In this paper, based on the theory of adaptive dynamics (Dieckmann and Law 1996; Geritz et al. 1998; Metz et al. 1992; Kisdi 2020; Hoyle et al. 2008; Geritz et al. 2007; Diekmann et al. 2005), we study the coevolutionary dynamics of host resistance and pathogen virulence. Based on an susceptible-infected population dynamics with density-dependent mortality, we will rigorously analyze the evolutionary invasion process of host resistance and pathogen virulence. In particular, strict mathematical proofs and numerical simulations are provided for each step of evolutionary invasion analysis. To sum up, the purpose of this study is aim to examine the following four questions of host-pathogen coevolution. Firstly, under what conditions will host and pathogen evolve to a continuously stable strategy? Secondly, under what conditions will the phenotypic traits of host and pathogen evolve periodically? Thirdly, under what conditions will the phenotypic traits of host or pathogen give rise to evolutionary branching? Fourthly, after evolutionary branching, can the host and pathogen with different phenotypic traits coexist stably for a long time?

The rest of the paper is organized as follows. In the next section, we first develop a susceptible-infected population model with density-dependent mortality. Then we use the approach of adaptive dynamics to analyze the evolutionary invasion process of host resistance and pathogen virulence and propose a coevolutionary dynamic model of host-pathogen interaction. In Sect. 3, the convergence stability, evolutionary stability and cyclical evolution are studied. After exploring the evolutionary branching of host resistance and pathogen virulence, the coevolutionary dynamics of two hosts and one pathogen and the coevolutionary dynamics of two pathogens and one host are given. In addition to the mathematical proof, the numerical simulation results are also given in Sect. 3 to illustrate the feasibility of our main results. A brief discussion is given at the end of this paper.

Evolutionary invasion analysis

In this section, we use the approach of adaptive dynamics to analyze the evolutionary invasion process of host and pathogen and establish a coevolutionary dynamic model of host-pathogen interaction. We first develop a susceptible-infected population model. Then based on this population dynamics, we will derive the invasion fitness for the mutant host and mutant pathogen that we use to investigate the coevolutionary dynamics and put forward to the proposition that successful invasion will lead to trait substitution. Finally, a coevolutionary dynamic model of host-pathogen interaction is constructed.

Susceptible-infected population dynamics

We assume that there is only one monomorphic host and one monomorphic pathogen in the initial community. Then we consider a pathogen (or a microparasite) that can cause a direct transmission with no latent period and the incidence of host infection is bilinear. In addition, when the host population density becomes larger, all hosts will be subject to intraspecific competition, the host mortality will increase with host density (Andreasen and Pugliese 1995), so we assume that both the susceptible host and infected host are subject to density-dependent mortality m(N), which is given by a linear function of the total population density , i.e.,where is the host natural death rate, is a constant representing the strength of density-dependence. In other words, the host population grows according to the logistic model. We also assume that infected individuals neither reproduce nor recover, which is a reasonable assumption for most invertebrate hosts with pathogen. In this way we focus on one possible route to disease resistance, that is, avoidance of infection by reducing the susceptibility to the disease. Therefore, the population dynamics of host-pathogen interaction is given bywhereand S and I denote the population density of susceptible host and infected host, respectively. is the host death rate due to the pathogen infection (i.e., the virulence), b is the birth rate of susceptible host, and is the transmission rate. All of the parameters in the model (1) are positive.

In the real world, the interaction between hosts and pathogens is commonly a coevolutionary process, the long-term behavior of host-pathogen interaction is likely to depend on the interplay of both species’ evolutionary characteristics (Best et al. 2009; Boots et al. 2009), so we assume that both the host’s resistance and the pathogen’s virulence can adaptively evolve. Host resistance and pathogen virulence are represented by a single phenotypic trait and , respectively. Specifically, on the one hand, the ubiquity of infectious pathogens in nature and the damage that they cause to their hosts has led to the evolution of a diverse range of host defence mechanisms, from simple mechanical barriers through to complex immune systems (Boots et al. 2009), so an increase in host resistance will lead to a decrease in transmission efficacy. On the other hand, although virulence is not an adaptation of the pathogen per se, it is generally believed to be an inevitable by-product of a pathogen’s need to propagate and transmit to new hosts, so an increase in virulence will parallel an increase in transmission efficacy (Anderson and May 1982). Hence, we assume that the transmission rate is an asymmetric function of trait , namely, the transmission rate decreases with the increase of the trait . To be specific, when the host resistance trait is much greater than the pathogen virulence trait , the transmission rate becomes very small. Conversely, when the pathogen virulence trait is far greater than the host resistance trait , then the transmission rate becomes very large. In addition, trade-offs are fundamental in nature and reflect the fact that any gain in one life-history trait (e.g., host resistance) incurs a cost in another (e.g., birth rate). Such trade-off relationship is also well documented (Boots and Begon 1993; Kraaijeveld et al. 2001; Boots et al. 2009). Boots and Begon (1993) examined genotypic trade-offs with resistance to a virus in a lepidopteran host by a micro-evolutionary selection experiment and found that correlated with the increase in resistance were a lengthening of development time, a reduction in egg viability and an increase in pupal weight. The biological basis of the resistance cost may be a direct pleiotropic effect of the alleles for resistance (Lenski 1988; Simms 1992). Hence, we further assume that a mutation that not only reduces the chance that an individual becomes infected but also includes a cost such that birth rate is reduced, i.e., the birth rate b is a monotonically decreasing function (trade-off function) with respect to the resistance trait , that is . Furthermore, the evolution of virulence can be understood as a balance between the pathogen’s drive to increase infectivity and its need to maintain a long infectious period. Effective transmission requires a high pathogen concentration in the host tissue while, on the other hand, a high pathogen production leads to extensive damage of host cells and consequently to shorter infection period (Anderson and May 1982). Hence, the increase in pathogen virulence will lead to an increase of the disease-related death rate, so we assume that the death rate due to disease is a monotonically increasing function (trade-off function) with respect to the virulence trait , that is, . To sum up, the population dynamics of host-pathogen interaction is changed towhereAll of the parameters in the model (2) are positive. This model is most applicable to invertebrate diseases and plant diseases (Boots and Haraguchi 1999).

Biologically, if the curve of trade-off function is concave, it means that as the resistance trait increases, will decrease more and more slowly. In this case, we say that there is a decelerating cost for the host population. On the contrary, if the curve of trade-off function is convex, which means that as the resistance trait increases, will decrease faster and faster, then we say that there is an accelerating cost for the host population. In the same way, if the curve of trade-off function is concave, in this case, we say that there is an accelerating cost for the pathogen population. If the curve of trade-off function is convex, then we say that there is a decelerating cost for the pathogen population. In addition, if the asymmetric function of transmission rate is concave, since is a decreasing function of , in this case, we say that there is a decelerating benefit for the host population, and an accelerating benefit for the pathogen population, and if the asymmetric function of transmission rate is convex, then we say that there is an accelerating benefit for the host population and a decelerating benefit for the pathogen population (Boots et al. 2009; Hoyle et al. 2008; Zu et al. 2015, 2020).

For model (2), the condition for the initial spread of a pathogen arriving into a susceptible host population is thenthat is, the basic reproduction numberWhen the condition (3) is satisfied, setting the right-hand side of the model (2) to zero, we obtain an endemic equilibrium , whereBy using the approach of Lyapunov function and the LaSalle’s invariance principle, we obtain the following result about the globally asymptotical stability of .

Proposition 1

If condition (3) holds, then the endemic equilibrium of the model (2) is globally asymptotically stable in .

Proof

For simplicity, we use instead of . Consider the following Lyapunov functionIt is clear that and the equality holds only for in . Its time derivative along the solution of model (2) becomesIt can be seen that in and if and only if . The globally asymptotical stability of follows from Lyapunov and LaSalle’s invariance principle.

Next, based on the susceptible-infected population dynamics (2) and endemic equilibrium , we derive the invasion fitness and the coevolutionary dynamics of host-pathogen interaction.

Invasion fitness and trait substitution

Based on the method of adaptive dynamics, we obtain the invasion fitness for mutant susceptible hostIf , then the population density of mutant host will increase (a detailed derivation is provided in Appendix A). In this case, we can say that the mutant host can invade. In addition, by using the approach of Lyapunov function, we find that the successful invasion of the mutant host will lead to a trait substitution (Cantrell et al. 2017; Dercole and Rinaldi 2008; Geritz et al. 2002; Geritz 2005; Meszéna et al. 2005) (a detailed proof is provided in Appendix B).

Analogously, the invasion fitness for mutant pathogen is given byIf , then the population density of mutant pathogen can invade (a detailed derivation is provided in Appendix A). By using the approach of Lyapunov function, we can see that the successful invasion of a mutant pathogen will result in a trait substitution (a detailed proof is provided in Appendix B).

Coevolutionary dynamics of host-pathogen interactions

Through successive invasion and trait substitution, the resistance of susceptible host and the virulence of pathogen will evolve step by step. The direction of such an evolutionary change is determined by the sign of selection gradient and , which is given bywhereWhen the mutation is small and rare, the coevolutionary dynamics of traits and is given bywhere time spans the evolutionary timescale, and are the selection gradients described as in (7), respectively, and represent the evolution speed of traits and , respectively, which is given bywhere and are the probability that a birth event in the host and pathogen is a mutant, respectively. and are the variance of phenotypic effect of the mutant host and mutant pathogen, respectively. Model (8) is an approximate equation of the coevolutionary process, which describes how the expected value of phenotypic traits will change.

Coevolutionary outcomes

Continuously stable strategy

When there is a pair of traits to satisfy the following conditionWe call for an evolutionary singular strategy (Geritz et al. 1998; Kisdi 1999). Whether the coevolutionary process of host-pathogen interaction stops at the singular strategy is determined by its convergence stability and evolutionary stability.

Convergence stability

We use the linear approximation method to estimate the convergence stability of the evolutionary singular strategy in model (8). The Jacobian matrix at the singular strategy is given byIf the determinant of this Jacobian matrix is positive , and the trace is negative , then the evolutionary singular strategy is locally convergence stable (Zu et al. 2014).

Evolutionary stability

The evolutionary singular strategy is evolutionarily stable means that it cannot be invaded by any nearby strategy, and it can be estimated by calculating the second derivative of the invasion fitness function of host and pathogen with respect to the mutant trait. By direct calculation, we obtain if the following condition is satisfiedwhereThen the evolutionary singular strategy is evolutionarily stable (Geritz et al. 1998; Zu et al. 2016).

If the evolutionary singular strategy is both convergence stable and evolutionarily stable, then the evolutionary singular strategy is a continuously stable strategy (CSS). Based on the above analysis, we obtain the following conclusion.

Theorem 1

Assuming condition (3) holds, for the evolutionary singular strategy of model (8), if , and condition (11) is satisfied, then the evolutionary singular strategy is a continuously stable strategy.

From the Jacobian matrix and condition (11), we can see that whether the singular strategy is a CSS depends on the curvature of the trade-off function at and the relative strength of the asymmetric interaction. At the same time, the population density at the endemic equilibrium of the susceptible host and the infected host also plays a key role. For the continuously stable strategy , the host and pathogen can stably coexist on a long-term evolutionary timescale. So a continuously stable strategy represents the final outcome of the coevolutionary process. In this case, the finally evolutionary outcome contains a monomorphic host and a monomorphic pathogen.

In order to give an example of numerical simulation to illustrate the result of Theorem 1, we take the following asymmetric transmission rate functionWhich is a monotonically decreasing function with respect to , where is the maximum transmission rate of pathogen, measures the strength of asymmetric interaction, the larger the , the stronger the asymmetric interaction, and adjusts the concavity and convexity. Particularly, when , is a turning point. When , if , then the transmission rate function is concave. That is, there is a decelerating benefit for the host population and an accelerating benefit for the pathogen population. If , then the transmission rate function is convex, that is, there is an accelerating benefit for the host population and a decelerating benefit for the pathogen population (see Fig. 1a). In general, this function can be suitable for a variety of asymmetric interaction (Kisdi 1999; Zu et al. 2016).

Fig. 1

Curves of transmission rate function and two trade-off functions. a Asymmetric curve of transmission rate function as given by (12), where . b The curve of birth rate function as given by (13), where . c The curve of disease-related mortality function as given by (14), where

In addition, based on the basic assumption, the birth rate b of the susceptible host population is a monotonically decreasing function with respect to the resistance trait , so we consider the birth rate function in the form ofWhen , the trade-off curve is globally convex, which means that the host population has an accelerating cost (see Fig. 1b). Similarly, the function form of the disease-related mortality is given byIt can be seen from Fig. 1c that the mortality function due to disease is a sigmoid curve with saturation, which is more in line with the actual biological significance, and the convexity is determined by the shape parameter , the larger the , the wider the convex area of the curve.

In order to illustrate the finally evolutionary outcome of host-pathogen coevolution, we take as bifurcation parameters in the following numerical simulation, and fix other parameters: . Below we give a numerical simulation example to illustrate the continuously stable strategy of host-pathogen coevolution. When , and , there is an evolutionary singular strategy , which is both convergence stable and evolutionarily stable, so represents the finally evolutionary outcome of host-pathogen coevolution. In particular, when , from Fig. 2a, it can be seen that is convergence stable. From Figs. 2b and 2c, it can be seen that the invasion fitness at reaches to a maximum value, so is evolutionarily stable. Therefore, is a continuously stable strategy. In this case, the finally evolutionary outcome contains a monomorphic host and a monomorphic pathogen, which can continuously stably coexist for a long time (see Fig. 2d). From Figs. 2e and 2f, we can see that the equilibrium population density of susceptible host and infected host finally reaches to a stable state, and the equilibrium population density of infected host is higher than that of the susceptible host.

Fig. 2

Continuously stable strategy of host-pathogen coevolution. a Trait evolution phase plot. The vector fields obtained from Model (8) indicate directions of coevolution of traits and . The black curve and red curve indicate respectively isoclines of traits and . The solid curves indicate the evolutionary singular strategy which is evolutionarily stable, while the dashed curve indicates the evolutionary singular strategy, which is not evolutionarily stable. The grey region is a feasible phenotypic trait space in which the host-pathogen coevolution can occur. b Invasion fitness landscape of the mutant host. c Invasion fitness landscape of the mutant pathogen. d Time series curves of phenotypic traits obtained through simulation of Model (8) with initial values are . e Equilibrium population density of the susceptible host when the traits and evolve. f Equilibrium population density of the infected host when the traits and evolve. Parameter values: . Other parameter values are the same as in Fig. 1

Evolutionary cycling

Based on the method of linear approximation, we can see that the convergence stability of evolutionary singular strategy is estimated by the determinant and trace of Jacobian matrix . If the strength of the trade-off function of disease-related mortality becomes stronger, then the evolutionary singular strategy may become unstable. In this case, the model (8) may admit Hopf bifurcation. That is, the phenotypic trait of host and pathogen may evolve to a stable limit cycle. So the evolutionary cycling is a possible outcome of host-pathogen coevolution. Due to the complex nonlinearity, the result of evolutionary cycling is illustrated by a numerical simulation example.

As an example, when , and , the evolutionary singular strategy becomes unstable. In this case, we find that the model (8) admits Hopf bifurcation, and the trait of the host and pathogen will converge to a limit cycle. The magnitude of the limit cycle increases with the increase of and finally disappears when reaches the critical value of 8.1. Therefore, we can see that the evolutionary cycling may occur under asymmetric interaction between the host and pathogen. In particular, when , that is, the strength of the trade-off function becomes stronger, from Fig. 3, we can see that the host resistance trait and the pathogen virulence trait will evolve periodically. In this case, the equilibrium population density of host and pathogen will also change periodically.

Fig. 3

Evolutionary cycling. a Trait evolution phase plot. The vector fields obtained from Model (8) indicate directions of coevolution of traits and . The black curve and red curve indicate respectively isoclines of traits and . The solid curves indicate evolutionary singular strategy which is evolutionarily stable, while the dashed curve indicates evolutionary singular strategy, which is not evolutionarily stable. The grey region is a feasible phenotypic trait space in which the host-pathogen coevolution can occur. b Time series curves of phenotypic traits obtained through simulation of Model (8) with initial condition . Host and pathogen evolve to a stable limit cycle. c Equilibrium population density of the host when the trait and evolve. d Equilibrium population density of the pathogen when the trait and evolve. Parameter values: . Other parameter values are the same as in Fig. 1

Evolutionary branching

Evolutionary branching is an evolutionary process during which directional selection drives a monomorphic population to an evolutionary singular strategy where ecological interactions induce disruptive selection and subsequently splits up into two coexisting phenotypic clusters (Doebeli and Dieckmann 2000). In one-dimensional adaptive dynamics, a fitness minimum at an attracting singular strategy is sufficient for evolutionary branching (Geritz et al. 1998). However, in higher-dimensional adaptive dynamics, this need not be true anymore. In higher-dimensional adaptive dynamics, an ‘evolutionarily singular coalition’ that is convergence stable but for which at least one strategy lacks evolutionary stability and allows for mutual invasibility nearby will lead to evolutionary branching (Dieckmann and Doebeli 1999; Doebeli and Dieckmann 2000; Kisdi 1999). Therefore, we obtain the following result on the evolutionary branching of the host and pathogen.

If the evolutionary singular strategy is convergence stable, and the evolutionary singular strategy of pathogen virulence is evolutionarily stable, but the evolutionary singular strategy of host resistance is not evolutionarily stable and allows for mutual invasibility nearby, which is given bythen the evolutionary branching of the host resistance will occur.

Similarly, we obtain if the evolutionary singular strategy is convergence stable, and the evolutionary singular strategy of host resistance is evolutionarily stable, but the evolutionary singular strategy of pathogen virulence is not evolutionarily stable and allows mutual invasibility nearby, which is given bythen the evolutionary branching of the pathogen virulence will occur.

To sum up, we obtain the following result about the evolutionary branching.

Theorem 2

Assuming condition (3) holds, for the evolutionary singular strategy of model (8),

(I) if , , and conditions (15), (16) and the second condition of equation (11) are satisfied, then the evolutionary branching of host resistance will occur;

(II) if , , and conditions (17), (18) and the first condition of Eq. (11) are satisfied, then the evolutionary branching of pathogen virulence will occur.

Based on the above analysis, we can see that whether the evolutionary branching of the host or pathogen occurs depends on the shape and relative strength of the transmission rate function and the strength of the trade-off function. From the condition (15) and (17), it can be seen that the equilibrium population density of the infected host and susceptible host also has a certain influence on whether the evolutionary branching occurs or not. In particular, from condition (15) it can be seen that if the birth rate function of the susceptible host is weakly convex at , but the transmission rate function is relatively strongly convex at , that is, the host has a weakly accelerating cost in term of birth rate, but at the same time has a relatively strongly accelerating benefit in term of transmission rate, then the evolutionary branching of host resistance may occur. In addition, from condition (17) we can see that if the disease-related mortality function is relatively strongly convex at , but the transmission rate function is weakly convex at , in other words, the pathogen has a weakly decelerating benefit in term of transmission rate, but at the same time has a relatively strongly decelerating cost in term of disease-related mortality, then the evolutionary branching of pathogen virulence may occur.

If the evolutionary branching of host resistance and pathogen virulence occurs, the host and pathogen will firstly evolve toward an evolutionary singular strategy , in the vicinity of the evolutionary singular strategy the host and pathogen will branch into two different types. After the evolutionary branching, we will further study the coevolutionary dynamics of two hosts and one pathogen or one host and two pathogens and explore the finally evolutionary outcomes.

Coevolutionary dynamics of two hosts and one pathogen

After the evolutionary branching of host resistance first occurs, we assume that there are two different types of hosts with resistance traits and , then the population dynamics of two-hosts-one-pathogen interaction is given bywhere .

Setting the right-hand side of (19) to 0, and letWhen the following condition is satisfiedwe obtain the endemic equilibrium of model (19), whereand .

By using the method of Lyapunov function, we will show that the endemic equilibrium is globally asymptotically stable.

Proposition 3

If condition (20) holds, then the endemic equilibrium , of model (19) is globally asymptotically stable in .

For simplicity, we use instead of . Consider the following Lyapunov functionIt is clear that and the equality holds only for in . Its time derivative along the solution of model (19) becomesIt can be seen that in and if and only if . The globally asymptotical stability of follows from Lyapunov-LaSalle’s invariance principle.

Next, based on model (19) and the endemic equilibrium, we will derive the coevolutionary dynamics of two hosts and one pathogen. Due to the rarity of mutation, we assume that there is either a mutant host arising from susceptible host or a mutant host arising from susceptible host , but not both at a time. By using the same derivation as before, when a mutant host with a different trait enters into the resident community with a low density, the invasion fitness for the mutant host is then given bywhere is the trade-off function of birth rate and are respectively the equilibrium population density of susceptible host , susceptible host and infected host I, which are described as in (21).

Analogously, the invasion fitness of mutant pathogen with different trait is given byTherefore, the evolutionary direction of host resistance and pathogen virulence is determined by the sign of selection gradients , and , which are given by

Furthermore, if the mutation is small and rare, then the coevolutionary dynamics of phenotypic traits is given bywhere and are the selection gradients described in (24), and are the probability that a birth event in the host and pathogen is a mutant, respectively. and are the variance of the phenotypic effect of mutant host and mutant pathogen, respectively. Model (25) is an approximate equation of coevolutionary process, which tells us how the expected value of traits will change.

From (24), we can get the evolutionary singular strategy . The convergence stability of the evolutionary singular strategy can be studied by using numerical simulation. The evolutionary stability depends on the following conditionTherefore, based on the above analysis, we obtain the following conclusion.

Theorem 3

Assuming condition (20) holds, for the evolutionary singular strategy of model (25), if it is convergence stable and satisfies the condition (26), then is a continuously stable strategy.

From the condition (25) and (26), it can be seen that whether the singular strategy is continuously stable or not depends on the strength and shape of the asymmetric host-pathogen interaction and the curvature of the trade-off function. It also depends on the equilibrium population density of susceptible host , susceptible host and the infected host I. Combining the previous analysis of the evolutionary branching of host resistance, we can see that if the host has a weakly accelerating cost in term of birth rate, but at the same time has a relatively strongly accelerating benefit in term of transmission rate, then the evolutionary branching of host resistance may occur. After the evolutionary branching, if condition (26) is satisfied, two hosts and one pathogen may converge to an evolutionary stable equilibrium. If the evolutionary singular strategy is both convergence stable and evolutionarily stable, then the singular strategy is the finally evolutionary outcome. In this case, two hosts and one pathogen can coexist stably for a long time. If the evolutionary singular strategy is convergence stable, but one singular strategy of the two hosts is not evolutionarily stable and allows for mutual invasion, then the host population will undergo further evolutionary branching. In this case, by using the same method as above, we can continue to study the final outcome of the further coevolution of multiple hosts and one pathogen.

Below we give a numerical simulation example to illustrate the evolutionary branching of host resistance. When , and , that is, the asymmetric host-pathogen interaction becomes relatively strong, we find that there is an evolutionary singularity strategy , which is convergence stable, but the host singulary strategy is not evolutionarily stable, so the evolutionary branching of the host may occur. Particularly, when , we can see that the evolutionary singular strategy . In this case, the host and pathogen will firstly evolve towards the singular strategy (0.320, 0.567). Near this evolutionary singular strategy, the asymmetric host-pathogen interaction becomes stronger, but there is a relatively strongly accelerating benefit, so the host population will branch into two different types (see Fig. 4). It can be seen from Fig. 4a that the evolutionary singular strategy is convergence stable. From Figs. 4b and 4c, it can be seen that the invasion fitness of the mutant host population minimizes at the resistance trait , so the resistance trait is not evolutionarily stable, and the mutant host and the resident host can invade each other, so the host resistance will undergo evolutionary branching. After the evolutionary branching of host resistance, we further study the coevolutionary dynamics of two hosts and one pathogen. From Figs. 4d and 4e, we find that the three populations finally converge to an evolutionary stable strategy , and they will be able to stably coexist for a long time. In addition, from Figs. 4f and 4g,we find that the equilibrium population density of the susceptible host with a higher resistance is relatively larger.

Fig. 4

Evolutionary branching of host resistance. a Trait evolution phase plot. The vector fields obtained from model (8) indicate directions of coevolution of traits and . The black curve and red curve indicate respectively isoclines of traits and . The solid curve indicates the evolutionary singular strategy, which is evolutionarily stable, while the dashed curve indicates the evolutionary singular strategy, which is not evolutionarily stable. The grey region is a feasible phenotypic trait space in which host-pathogen coevolution can occur. b Pairwise invasibility plot for fixed pathogen strategy . c Mutual invasibility plot for fixed pathogen strategy . d Invasion fitness landscape when . e Simulated evolutionary tree obtained through simulation of model (8) with initial condition and model (25). f Equilibrium population density of susceptible host when the traits and evolve. g Equilibrium population density of infected host when the traits and evolve. Parameter values: . Other parameter values are the same as in Fig. 1

However, if , but , that is to say, the density-dependent natural mortality becomes relatively strong, then the host regains the evolutionary stability, and no evolutionary branching occurs. Instead, there is a continuously stable strategy . From Fig. 5a, it can be seen that is convergence stable. From Figs. 5b and 5c, it can be seen that the invasion fitness at reaches the maximum value, so is evolutionarily stable. Therefore, is a continuously stable strategy. The host and the pathogen can stably coexist for a long time. Comparing Fig. 4 with Fig. 5, we can see that only changing the value of results in two completely different evolutionary phenomena. The intensity of the density-dependent natural mortality greatly influences the coevolutionary outcome of the host and pathogen. In this case, the finally evolutionary outcome includes a monomorphic host and a monomorphic pathogen (see Fig. 5).

Fig. 5

Continuously stable strategy of host and pathogen when density-dependent natural mortality increases. a Trait evolution phase plot. The vector fields obtained from model (8) indicate directions of coevolution of traits and . The black curve and red curve indicate respectively isoclines of traits and . The solid curve indicates the evolutionary singular strategy, which is evolutionarily stable. The grey region is a feasible phenotypic trait space in which the host-pathogen coevolution can occur. b Invasion fitness landscape of the mutant host. c Invasion fitness landscape of the mutant pathogen. d Time series curves of phenotypic traits and obtained through simulation of model (8) with initial value . e Equilibrium population density of susceptible host when the traits and evolve. f Equilibrium population density of infected host when the traits and evolve. Parameter values: . Other parameter values are the same as in Fig. 1

Coevolutionary dynamics of two pathogens and one host

Analogously, after the evolutionary branching of pathogen virulence first occurs, we assume that there are two different types of pathogens with virulence traits and , then the population dynamics of the two pathogens and one host is given bywhere .

Setting the right-hand sides of (27) to 0, and letWhen the following condition is satisfiedwe obtain the endemic equilibrium of model (27), whereand .

Similar to the proof of Proposition 3, we can prove that the endemic equilibrium is globally asymptotically stable.

When a mutant host with a different trait enters into the resident community with a low density, the invasion fitness for the mutant host is given byAnalogously, the invasion fitness for the mutant pathogen with different trait is given byThus, the coevolutionary dynamics of traits , and is given bywhere and are the probability that a birth event in the host and pathogen is a mutant, respectively. and are the variance of the phenotypic effect of mutant host and mutant pathogen, respectively. Model (32) is an approximate equation of coevolutionary process, which describes how the expected value of traits will change.

From (33), we can get the evolutionary singular strategy . Because of the very complex nonlinearity, we use the method of numerical simulation to estimate the convergence stability of singular strategy. The evolutionary stability depends on the following conditionBased on the above analysis, we obtain the following conclusion.

Theorem 4

Assuming condition (28) holds, for the evolutionary singular strategy of model (32), if it is convergence stable and satisfies the condition (34), then is a continuously stable strategy.

From the conditions (32) and (34), we can see that whether the evolutionary singular strategy is continuously stable not only depends on the strength and shape of asymmetric host-pathogen interaction, but also depends on the equilibrium population density of the susceptible host S and the infected hosts and . Particularly, combining with the previous analysis of the evolutionary branching of pathogen virulence, we can see that if the pathogen has a weakly decelerating benefit in term of transmission rate, but at the same time has a relatively strongly decelerating cost in term of disease-related mortality, then the evolutionary branching of pathogen virulence may occur. After the evolutionary branching occurs, if condition (34) is satisfied, then the two pathogens and one host may converge to an evolutionary stable strategy. If the singular strategy is both convergence stable and evolutionarily stable, then the singular strategy is the finally evolutionary outcome. In this case, the finally evolutionary outcome includes a dimorphic pathogen and a monomorphic host. Furthermore, if the singular strategy is convergence stable, but one singular strategy of the two pathogens is not evolutionarily stable and allows for mutual invasion, then the pathogen will further undergo evolutionary branching. In this case, we can use the same method as above to study the coevolutionary dynamics of multiple pathogens and one host.

Below we give a specific numerical simulation example to illustrate the evolutionary branching of pathogen virulence. When , and , that is, the strength of disease-related mortality becomes relatively strong, we find that there is an evolutionary singularity strategy , which is convergence stable, but the evolutionary singularity strategy of pathogen is not evolutionarily stable, so the evolutionary branching of pathogen virulence may occur. Particularly, when , we can see that there is an evolutionary singular strategy . It can be seen from Fig. 6a that the evolutionary singular strategy is convergence stable. From Figs. 6b and 6c, it can be seen that the virulence trait is convergence stable but not evolutionarily stable, and the mutant pathogen and the resident pathogen can invade each other, so the pathogen virulence will undergo evolutionary branching. In this case, the host and the pathogen will first evolve towards the singular strategy (0.396, 0.549). Near this evolutionary singular strategy, the strength of disease-related mortality becomes stronger, but the asymmetric interaction between the host and the pathogen becomes not very strong, that is, there is a relatively strongly decelerating cost in term of disease-related mortality for the pathogen population, so the pathogen virulence will branch into two different types (see Fig. 6). After the evolutionary branching of pathogen virulence, from Figs. 6d and 6e, we find that the three populations finally converge to an evolutionary stable strategy , and they can continuously stably coexist for a long time. From Figs. 6f and 6g, it can be seen that the equilibrium population density of the two pathogens and the equilibrium population density of the one host finally reached to a stable state.

Fig. 6

Evolutionary branching of pathogen virulence. a Trait evolution phase plot. The vector fields obtained from model (8) indicate directions of coevolution of traits and . The black curve and red curve indicate respectively isoclines of traits and . The solid curve indicates the evolutionary singular strategy which is evolutionarily stable, while the dashed curve indicates the evolutionary singular strategy which is not evolutionarily stable. The grey region is a feasible phenotypic trait space, in which the coevolution of host and pathogen can occur. b Pairwise invasibility plot for fixed host strategy . c Mutual invasibility plot for fixed host strategy . d Invasion fitness landscape when . e Simulated evolutionary tree obtained through simulation of model (8) with initial condition and model (32). f Equilibrium population density of the infected host with a pathogen when the traits and evolve. g Equilibrium population density of the susceptible host when the traits and evolve. Parameter values: . Other parameter values are the same as in Fig. 1

Combined with Figs. 2– 6, we find that the finally coevolutionary outcome of host-pathogen interaction not only depends on the strength of the asymmetric interaction, but also depends on the strength of density-dependent natural mortality and the intensity of the two trade-off functions. The finally coevolutionary outcome of host-pathogen interaction may include continuously stable strategy, evolutionary cycling, evolutionary branching of host resistance, evolutionary branching of pathogen virulence, and continuously stable coexistence of two-hosts-one-pathogen and one-host-two-pathogens.

Discussion

Since that most natural interaction system will be coevolutionary, therefore, it is essential to develop a basic framework of coevolutionary dynamics to completely understand the evolutionary mechanism and outcome. In this paper, we study the coevolutionary dynamics of host-pathogen interaction in detail. Specifically, we examine four questions of host-pathogen coevolution. Firstly, under what conditions will host and pathogen evolve to a continuously stable strategy? Secondly, under what conditions will host and pathogen phenotypic traits evolve periodically? Thirdly, under what conditions will host and pathogen phenotypic traits give rise to evolutionary branching? Fourthly, after evolutionary branching, can the hosts and pathogens with different phenotypic traits coexist stably for a long time? Based on the theory of adaptive dynamics and numerical simulation, we find that the evolutionary diversity of host resistance and pathogen virulence is driven by multiple factors, not only depends on the shape and strength of the host-pathogen asymmetric interaction, but also depends on the intensity of the two trade-off functions and the equilibrium population density of susceptible host and infected host. In general, this research reaches important conclusions in the following three aspects. Firstly, it can be seen that if the strengths of host-pathogen asymmetric interaction and disease-related mortality are relatively weak, then the host resistance and pathogen virulence will evolve to a continuously stable strategy. However, if the strength of host-pathogen asymmetric interaction and disease-related mortality becomes stronger, then the host resistance and pathogen virulence will evolve periodically. Secondly, we find that if the birth rate function of susceptible host is weakly convex at the evolutionary singular strategy, but the transmission rate function is relatively strongly convex at the evolutionary singular strategy, that is, the host has a weakly accelerating cost in term of birth rate, but at the same time has a relatively strongly accelerating benefit in term of transmission rate, then the evolutionary branching of host resistance may occur firstly. However, if the disease-related mortality function is relatively strongly convex at the evolutionary singular strategy, but the transmission rate function is weakly convex at the evolutionary singular strategy, in other words, the pathogen has a weakly decelerating benefit in term of transmission rate, but at the same time has a relatively strongly decelerating cost in term of disease-related mortality, then the evolutionary branching of pathogen virulence may occur firstly. Thirdly, after the evolutionary branching of host resistance and pathogen virulence, we further study the coevolutionary dynamics of two-hosts-one-pathogen interaction and one-host-two-pathogens interaction. We find that if the evolutionary branching of host resistance arises firstly, then the finally evolutionary outcome may contain a dimorphic host and a monomorphic pathogen population. If the evolutionary branching of pathogen virulence arises firstly, then the finally evolutionary outcome may contain a monomorphic host and a dimorphic pathogen population, which can stably coexist for a long time. These results are essential for our more general understanding of host-pathogen interaction and the coevolutionary mechanism.

Compared with the previous studies (Boots et al. 2009; Best et al. 2010, 2011; Lopez Pascua et al. 2014; Boots et al. 2014; McLeod and Day 2015; Hesse et al. 2015; Best 2018; Alizon 2021), the novelty of our work is reflected in the following several aspects: Firstly, in terms of models, previous studies only considered the density-dependent reproduction (Best et al. 2011; Boots et al. 2014), see also (Kada and Lion 2015). Based on a population model with density-dependent reproduction, Best et al. (2011) studied the impact of different degrees of spatial structure on the evolution of host resistance, and on the coevolution of host resistance and parasite virulence. They found that for the host, local reproduction will always select for higher resistance, and globalized interactions will select for diseases characterized by low host defenses, high disease transmission and high parasite virulence. Boots et al. (2014) examined the effect of epidemiological feedbacks and characteristics of interaction between host and parasite on the coevolution of host-parasite diversity. They showed that epidemiological feedbacks may generate diversity when host resistance and parasite infectivity have costs. For trait polymorphism, both specificity of infection between host and parasite, and incompatibility between particular strain and type are required. Kada and Lion (2015) showed that the result of coevolution depend critically on the specific ecological and biological ingredients of the system. In contrast with population models of density-dependent reproduction, our model assumes that the natural mortalities of both the susceptible host and infected host are density-dependent, which is more reasonable in the real world. In this case, we can also construct a Lyapunov function and prove the globally asymptotic stability of the endemic equilibrium and the result of ‘invasion implies trait substitution’. This proof is necessary for the rigorously evolutionary invasion analysis. Especially, we find that if the density-dependent natural mortality is relatively weak, then the evolutionary branching of host resistance can occur. However, if the density-dependent natural mortality becomes relatively strong, then the host regains the evolutionary stability, and no evolutionary branching occurs. Secondly, the two trade-off functions and the asymmetric transmission rate function take a more flexible form, which are more consistent with empirical evidence and suitable for a wider range of asymmetric host-pathogen interaction. Thirdly, we obtain the the ecological and evolutionary conditions for the evolutionary diversity of host resistance and pathogen virulence, and explore the dynamic evolutionary branching process by numerical simulation. After evolutionary branching, we further study the finally evolutionary outcome of host resistance and pathogen virulence. Fourthly, we perform the numerical bifurcation analysis on coevolutionary dynamics and discuss the effect of demographic parameter on the evolutionary behavior of host resistance and pathogen virulence, which is helpful for us to develop the bifurcation theory of evolutionary dynamics. Especially, the framework and method of this study also provide a foundation and technique for studying the evolutionary diversity of other infectious diseases.

However, considering the sensitivity of the interaction between the host and pathogen, we need to say (1) based on more complex epidemiological models, such as SIS or SIR models, studying the coevolutionary dynamics of host resistance and pathogen virulence will have more practical significance. (2) This paper considers that the host defense mechanism is to avoid infection, but the host defense mechanism is diverse, including not only avoiding infection but also recovering faster after infection or surviving longer once infected, so exploring the evolution of host resistance under different defense mechanisms will be a problem for our further research. (3) When performing numerical simulation on the result of coevolution, the form of each trade-off function is selected based on theoretical analysis and experience, and actual data fitting is not performed. Therefore, in future work, we will further collect actual data (such as, the mutation and variant data of SARS-CoV-2) and fit the most accurate trade-off function. (4) The coevolutionary process of the host and pathogen has not been fully discussed, and the host or the pathogen may have further evolutionary branching. In addition, we can explore the condition for the evolutionary extinction of one of the host and pathogen, which is worthy of our further exploration. Further research on these problems will not only help us further comprehensively understand the mechanism of host-pathogen interaction, but also play an important role in clarifying the evolutionary mechanism of biodiversity formation.