![MathJax Logo](/templates/jsp/_style2/_tandf/pb2/images/math-jax.gif)
ABSTRACT
In this paper, we use a two-host one pathogen immuno-epidemiological model to argue that the principle for host evolution, when the host is subjected to a fatal disease, is minimization of the case fatality proportion . This principle is valid whether the disease is chronic or leads to recovery. In the case of continuum of hosts, stratified by their immune response stimulation rate a, we suggest that
has a minimum because a trade-off exists between virulence to the host induced by the pathogen and virulence induced by the immune response. We find that the minimization of the case fatality proportion is an evolutionary stable strategy for the host.
1. Introduction
Pathogen evolution, and specifically the evolution of virulence, have been studied extensively through mathematical models with the early models assuming some kind of trade-off between virulence and transmission [Citation7,Citation11,Citation13,Citation20,Citation21,Citation25]. A competitive exclusion principle, suggesting that pathogens evolve to maximize their epidemiological reproduction number, was first rigorously established by Bremermann and Thieme [Citation12]. This principle lead the way into the following studies on pathogen evolution and to the fundamental idea that pathogens evolve to some optimal virulence, that maximizes the reproduction number.
Studies in evolutionary genetics show that hosts, including humans, in turn, evolve under the selective pressure exhorted by the pathogen [Citation6]. Mathematical models have investigated the evolution of host resistance [Citation8,Citation9,Citation16,Citation24]. Bowers [Citation9] derived a principle for host evolution, similar to the reproduction number. His principle states that the host evolves towards minimizing a dimensionless quantity called the basic depression number .
Anderson and May are the pioneers in the mathematical study of co-evolution of pathogens and hosts [Citation3]. Early models of co-evolution, reviewed in [Citation28], were mostly systems of ordinary differential equations (ODEs). In 2002, Gilchrist and Sasaki proposed a novel model of ‘nested’ dynamical immunological model into an epidemiological model [Citation17]. Since then nested immuno-epidemiological models have become a primary tool in the study of evolution of pathogens [Citation1,Citation2,Citation4,Citation5,Citation14,Citation15,Citation19]. Gilchrist and Sasaki, in fact, studied the co-evolution of the pathogen and the host. They used the lifespan of the host in the infectious class as a criterion for host evolution, assuming no natural mortality. Pugliese [Citation26] later showed analytically that, even in the presence of natural mortality, the host evolves towards maximizing its lifespan in the infectious class. This measure, however, is not adequate for immuno-epidemiological models that model diseases with recovery. Why would the host strive to extend its lifespan in the infectious class, if it has an option to recover?
The purpose of this note is to extend the criterion for host evolution to immuno-epidemiological models with recovery. We argue that the host evolves towards minimizing the case fatality proportion . The case fatality proportion in immuno-epidemiological models is given by
where
is the disease-induced death rate and
is the probability of survival in the infectious class. In the ODE case the case fatality proportion is given by
where γ is the recovery rate and
is the natural death rate. The case fatality proportion gives the fraction of the individuals who die from the disease. In epidemiology,
is called case fatality ratio (CFR). Our measure is in a sense a special case of the one derived by Bowers and assumes that all hosts have the same reproduction rate and the same natural death rate, that is, we assume that these host characteristics do not evolve under the pressure of the pathogen [Citation9]. We further consider the case fatality proportion
as a function of the immune response activation rate a and argue that the case fatality proportion has a minimum because a trade-off exists between the virulence generated by the pathogen and the virulence generated by the immune response. Using tools from adaptive dynamics, we show that minimizing the case fatality proportion is an evolutionary stable strategy (ESS) for the host.
This paper is organized as follows. In the next section, we introduce a two-host immuno-epidemiological model. We link the epidemiological parameters to the within-host dynamic variables [Citation17]. We derive the basic reproduction number. In Section 3 we compute the equilibria of the model and determine their stabilities. We show that host one equilibrium is stable if host one minimizes the case fatality proportion. If host one has the smallest case fatality proportion, the stability of host one equilibrium is established in several particular cases. In Section 4 we consider simulations with a continuum of hosts, stratified by their immune response stimulation rate and we show that
has a minimum. Section 5 summarizes our results.
2. A two-host immuno-epidemiological model
In this section we formulate a two-host immuno-epidemiological model which models directly transmitted diseases with recovery, such as Hepatitis C Virus, influenza and foot-and-mouth disease. Both hosts have the same recruitment rate and natural mortality but vary in susceptibility, recovery rate and disease-induced mortality (virulence). Both hosts experience different within-host dynamics, given by a standard immunological single-strain model.
2.1. The within-host model
Within a host, the time variable is denoted by τ and signifies time-since-infection. We model only the pathogen and the immune response, denoted, respectively, by
and
. The pathogen reproduces linearly at a rate
and is killed by the immune response with efficacy
. The pathogen stimulates the immune response at a rate
. The within-host model for host i takes the form [Citation17]:
(1)
(1)
This model is equipped with initial conditions and
.
All parameters and dependent variables of this within-host model and their definitions can be found in Table .
Table 1. List of parameters and variables.
This model has a very simple dynamics. Within a host, if the pathogen takes off, reaches a maximum and then declines to zero. The immune response increases and levels off. If
, the pathogen decreases to zero.
Dividing the two equations in (Equation1(1)
(1) ), we can obtain P as a function of B,
(2)
(2) Using the second equation in (Equation1
(1)
(1) ), we obtain the following differential equation for B:
(3)
(3) This is a separable equation but it is still hard to solve exactly because of the complexity of the right-hand side.
2.2. The between-host model
In this subsection we introduce a two-host epidemiological model. The two hosts are subjected to a unique pathogen. We denote by the number of susceptible individuals of host type j at time t. We structure the infected individuals by time-since-infection τ. Let
be the density of infected individuals of host type j. Individuals in the class
experience the same within-host dynamics given by model (Equation1
(1)
(1) ) for host type j. We assume the epidemiological dynamics is given by the following two-host model:
(4)
(4)
In model (Equation4(4)
(4) ),
is the birth/recruitment rate where
is the total population size of host j and
is the total population size of both hosts.
is a function of N such that
for all N and
exists. Furthermore,
. In addition,
is the natural death rate and
is the disease-induced death rate of host type j. We will assume the simplest form of
, that is,
where
represents host mortality due to the pathogen and
gives the additional host mortality due to the immune response. We note that other forms of disease-induced mortality rate are possible [Citation17]. The transmission coefficient
is also dependent on the within-host pathogen load. We may assume that
is Holling type II with respect to the pathogen load at a given time-since-infection τ. Hence,
where
is the viral load in host type j,
is the half-saturation constant and
is proportionality constant. The recovery rate is directly proportional to the immune response
and inversely related to the viral load. Hence,
where
is proportionality constant and
is a small number.
Model (Equation4(4)
(4) ) is equipped with the following initial conditions:
,
and
with j=1,2. All parameters are nonnegative and
,
.
The epidemiological reproduction number of host type j in system (Equation4(4)
(4) ) is given by the following expression:
(5)
(5) The reproduction number of host type j gives the number of secondary infections that one type j-infected individual will produce in an entirely susceptible type j population during its lifespan as infectious. This is the reproduction number of the single-host model and does not depend on the interaction between the hosts.
In the next section we compute explicit expressions for the equilibria and establish their local stability.
3. Equilibria and their local stability
System (Equation4(4)
(4) ) has several disease-free equilibria.
which always exists.
where
. This equilibrium exists if and only if
.
where
. This equilibrium exists if and only if
.
There is also one parameter family of coexistence equilibria
satisfying
Before we introduce the host type j equilibria, we define the probability for survival in the jth infected class,
If the basic reproduction number , for each j there is a corresponding single host infected equilibrium
given by
.
where
and
is the solution of the following equation:
(6)
(6) and
is given by the following constant:
This claim merits justification. We formulate the result in the following proposition:
Proposition 3.1.
Assume . Then there exists a unique host type j equilibrium
where
and
.
Proof.
We need to show that Equation (Equation6(6)
(6) ) has a unique solution satisfying
. If
, then the expression in the parenthesis is positive and, therefore,
. Hence the equilibrium
has only positive components for host type j. Given the properties of b, it is not hard to see that Equation (Equation6
(6)
(6) ) always has a unique positive solution.
To see that , we rewrite Equation (Equation6
(6)
(6) ) as
where we note that
. The left-hand side of this equation is a function of
, say
. Since
is an increasing function of
, if the above equation has a solution, it is unique. To see the above equation has a solution satisfying
, notice that
which follows from the fact that
and therefore,
. Let
be such that
. Then,
Hence, Equation (Equation6
(6)
(6) ) has a unique solution in the interval
. This concludes the proof.
Next, we investigate the local stability of equilibria. First, we consider the extinction equilibrium. In the case of the extinction equilibrium, we state the following result which is not hard to prove.
Proposition 3.2.
The extinction equilibrium is locally asymptotically stable, if
and unstable if
.
The stability of the disease-free boundary equilibria is given by the following proposition whose proof will also be omitted.
Proposition 3.3.
Equilibrium is neutrally stable
principle eigenvalue is zero
if
and unstable if
.
We state the stability of the one parameter family of coexistence equilibria in the following proposition. We first denote
Proposition 3.4.
Equilibrium is locally asymptotically stable if
and unstable if
.
Next we establish the stability of . Stability of
is similar. The main result, the stability of
, gives conditions for the outcome of the competition of hosts one and two. It is well known that the outcome of the competition of multiple strains where competitive exclusion is the only outcome, is governed by the reproduction number – the strain with the maximal reproduction number eliminates the rest [Citation12]. Here we establish that the competition between hosts, subject to the same pathogen, is governed by the case fatality proportion, defined as follows:
(7)
(7) We note that
. The case fatality proportion governs the outcome of the competition in both chronic diseases and diseases leading to recovery. Pugliese [Citation26] suggested a different measure that governs the outcome of host competition, the lifespan in the infected class, but this measure is valid only if
, that is, recovery is impossible.
Stability of will give the criterion for the outcome of host competition. However, stability
depends on two components (1) internal stability of the endemic equilibrium of a host one only system and (2) criterion for non-invasion of host two. In system (Equation4
(4)
(4) ) what is difficult to prove and may not hold for all parameter values is (1). In the theorem below, we first establish the criterion for non-invasion, assuming (1).
Theorem 3.5.
Assume the endemic equilibrium of host one only system is locally asymptotically stable. Then, the equilibrium is locally asymptotically stable if and only if
that is, the host with the smallest case fatality proportion outcompetes the rest.
Proof.
Denote . We linearize system (Equation4
(4)
(4) ) around equilibrium
by setting
,
,
,
,
,
,
and
. We look for exponential solutions. That leads to the following eigenvalue problem.
(8)
(8) We consider first the eigenvalues associated with competitor two. Solving the differential equation, we have
. Replacing
with
we obtain the following characteristic equation from the equation for
:
(9)
(9) Note that from the equilibrium equation for
, we have
Furthermore, from the equation for the equilibrium , we have
Hence the constant of the left-hand side of Equation (Equation9
(9)
(9) ) becomes
(10)
(10) Now, denote by
the left-hand side of Equation (Equation9
(9)
(9) ) and by
the right-hand side. If λ is a real variable, then
is a linear increasing function of λ and
is a decreasing function, approaching zero.
has a positive real solution if and only if
, that is if and only if
Thus, we conclude that the characteristic equation (Equation9
(9)
(9) ) has positive real root if and only if
and the equilibrium
is unstable in this case.
If the characteristic equation (Equation9
(9)
(9) ) does not have a positive real root. We show that it does not have complex roots with nonnegative real parts. Assume
and
. Then
At the same time
Since
, then
and
Thus the characteristic Equation (Equation9
(9)
(9) ) does not have roots with nonnegative real part. If
, stability of
depends on the eigenvalues of the system for
,
and
. If λ is distinct from the eigenvalues of Equation (Equation9
(9)
(9) ), then
and consequently
and
. The system for
,
and
is the system for the endemic equilibrium in the host one only model. Since we assume that equilibrium is locally asymptotically stable, that system have only eigenvalues with negative real part. Hence,
is locally asymptotically stable. That concludes the proof.
In the next theorem, we establish several conditions for stability of the endemic equilibrium of host one only system.
Theorem 3.6.
Assume . The endemic equilibrium of host one only system is locally asymptotically stable in one of the following cases.
that is
.
and
that is
.
and
. That is all linking parameters are set to constants.
Proof.
The linearized system for the endemic equilibrium of host one only takes the form:
(11)
(11)
From the equation for , we have
Then,
From the equation for
, we obtain
(12)
(12) where
. Replacing
in the equation for
, we obtain the following characteristic equation:
(13)
(13)
We consider three special cases.
Case 1 , that is
. Assume in addition that
. Then,
On the other hand, if
with
:
(14)
(14) That is, the left-hand side of the characteristic equation (Equation13
(13)
(13) ) is strictly bigger in absolute value than one, while the right-hand side is smaller or equal than one. Thus, these cannot be equal for complex λ with nonnegative real part. In this computation, we must recall that we have assumed
. Furthermore, we must mention that above sequence of inequalities is valid since
because
is a positive decreasing function of τ whose value at zero is positive.
Our inability to prove stability in the case may appear natural in the light of results of Thieme and Castillo-Chavez [Citation27]. However, nested immuno-epidemiological models tend to be a lot more stable than single scale age-since-infection structured models, as we have shown in [Citation23]. The component that can potentially destabilize the model here is the logistic recruitment function. If recruitment is constant, as in [Citation23,Citation27], then D=0 and the proof in Case 2 below works without the need to assume
.
Case 2 When and
, that is
. In this case, keeping in mind that
we can reduce Equation (Equation13
(13)
(13) ) to the following characteristic equation:
(15)
(15) Since
, we have
As before, it can be shown that this equation does not have roots with nonnegative real part since for complex λ with nonnegative real part in absolute value the left-hand side is bigger than one while the right-hand side is smaller than one.
Case 3 In this case we assume that ,
, and
. That is all linking parameters are set to constants and do not depend on the within-host model. This makes the epidemic model independent of the within-host model and turns it into an ODE. In this case, the characteristic equation (Equation15
(15)
(15) ) becomes
Collecting terms, this can be reduced to the following cubic equation:
where
(16)
(16)
Clearly, since
,
,
and
. It is not hard to check that
. Thus the Routh–Hurwitz conditions imply stability.
Remark 1.
In the case when ,
, and
,
,
, and
the condition
takes the form
(17)
(17)
We were not able to establish the local stability of the endemic equilibrium of host one in its general case. It is not hard to see that the characteristic equation (Equation15(15)
(15) ) does not have any nonnegative real roots. We believe, however, that this equation may have complex roots with positive real parts. All special cases that we considered in an attempt to find complex roots with positive real part resulted in special cases of the characteristic equation with only roots with negative real parts. Thus, we were not able to destabilize the system and obtain oscillations, although we could not rule those out.
4. Minimizing the case fatality proportion
Our results in the previous section suggest that on population level the host evolves by minimizing the case fatality proportion. This principle is valid both in diseases without recovery and in diseases with recovery. Gilchrist and Sasaki [Citation17] suggested that on within-host level the host evolves by changing its immune response activation rate a. Thus, instead of considering two hosts, we consider a continuum of hosts, whose case fatality proportion is a function of the immune response activation rate a. The host evolves towards minimizing
. Why should
have a minimum? We surmise that as a increases
decreases as the better immune response suppresses the pathogen. Thus, the case fatality proportion would decrease if the virulence were only generated by the pathogen. However, as a increases the immune response
increases and the case fatality proportion increases in the absence of pathogen-induced mortality (see Figure ). Thus, from the view point of the host, a trade-off exists between virulence induced by high pathogen load and virulence induced by high immune response. The host has to optimize its immune response so that it controls the pathogen but does not affect the host too much.
Table 2. List of Parameter Values.
Pathogens evolve very quickly and in the process they change their reproductive rate r. The host has to adapt to the evolving pathogens. How do the case fatality proportion and the optimal immune response change with the pathogen changing reproductive rate? Figure suggests that as r increases, the case fatality proportion also increases. This observation is somewhat intuitive as increasing r increases , which increases disease-induced mortality, and ultimately increases the case fatality proportion. We expect that this observation is parameter-dependent and different outcome may be possible.
Figure 2. The case fatality proportion as a function of a for three values of r. Black points denote the minimum. Parameters are given in Table .
![Figure 2. The case fatality proportion as a function of a for three values of r. Black points denote the minimum. Parameters are given in Table 2.](/cms/asset/42ec8517-0bed-43e8-a7fe-0a70a122b14f/tjbd_a_1161089_f0002_c.jpg)
Figure also suggests that as r increases, the optimal immune response rate also increases to compensate for the increased pathogen reproduction. This observation seems robust and was also made in the case when host lifespan in the infected class was considered as host optimization principle [Citation17,Citation26].
Just like maximizing the epidemiological reproduction for the pathogen is an ESS, minimizing the case fatality proportion for the host is also an ESS. ESS is a concept from adaptive dynamics which studies the evolution of the traits. Adaptive dynamics considers the long-term consequences of the potential invasion of a mutant trait, when the resident population, which adopts a resident trait, is at equilibrium. In our setting here the trait under evolution is the immune activation rate a. We will use the case fatality proportion as a proxy to the invasion fitness [Citation10]. In the case of pathogen evolution this role is played by the reproduction number. Actually, the case fatality proportion in a sense is ‘anti-fitness’. The fitness for the host will be . Nonetheless, we will work with the case fatality proportion. We know that the case fatality proportion is function of a:
. Traits for which
are called evolutionarily singular strategies. We know that the case fatality proportion can have a minimum, where the derivative vanishes. Hence, there is at least one evolutionary singular strategy. Evolutionary singular strategy can be ESS, branching point or a singular case. That can be determined graphically from a pairwise invasibility plot (PIP). PIP gives the region where the mutant using a mutant strategy
can invade the resident, playing resident strategy
. That region is obtained from plotting the region where the inequality
holds in the
plane (see Figure ). The evolutionarily singular strategies in a PIP are obtained from the intersection of the boundary of the invasion region with the main diagonal. Figure has one evolutionarily singular strategy. This corresponds to the minimum of
in Figure . The evolutionarily singular strategy is an ESS if the vertical line that passes through the singular strategy lies entirely in the region of non-invasion. That signifies that the evolutionarily singular strategy, once established, cannot be invaded by nearby mutants, that is it is an ESS. That is the case in Figure . Hence, minimizing the case fatality proportion is an ESS for the host.
Figure 3. Pairwise Invasibility Plot: The x-axis is the immune activation rate of the resident host. The y-axis is the immune activation rate of the mutant host. The green colour represents the region where .
![Figure 3. Pairwise Invasibility Plot: The x-axis is the immune activation rate of the resident host. The y-axis is the immune activation rate of the mutant host. The green colour represents the region where F(am)<F(ar).](/cms/asset/eb8d63b1-a465-46a6-b0d2-52c4f9e44199/tjbd_a_1161089_f0003_c.jpg)
We would like to mention that although we do recognize pathogens will generally evolve faster than hosts, so that the hosts have to adapt to evolving pathogens, this aspect is not considered in building the PIP, which is based on a fixed parameter value for the pathogen. The reason for that is that in this note we are only interested in the evolution of the host as a response of an ‘evolutionary fixed’ pathogen, and not in the host–pathogen co-evolution.
One interesting observation, suggested by a referee, is that in the current immuno-epidemiological model, the ‘fitness’ of an invader host does not depend on the resident host type. In such cases, it is clear that ESS will minimize the case fatality rate, and more complex evolutionary outcomes cannot occur. Indeed, the PIP is perfectly symmetrical (if , host 2 cannot invade host 1 while host 1 can invade host 2).
5. Discussion
In this paper we consider a two-host single pathogen immuno-epidemiological model. The two hosts are subjected to competitive exclusion. The main point of the paper is to show that the host that has the smallest population level case fatality proportion persists, while the other dies out. Without loss of generality, we assume that host one has a smaller case fatality proportion
, and we consider local stability of equilibrium
. We show that if
, host two cannot invade host one and the stability of
is further determined by the stability of the characteristic equation for host one in the absence of host two. We were able to establish stability of this equation in three special cases: (1) when recovery rate is zero and the derivative of the recruitment function at
is nonnegative and smaller that the natural death rate; (2) when the disease-induced death rate is zero; (3) when all epidemiological parameter functions are constants and if the derivative of the recruitment function at
is nonnegative, it is smaller that the natural death rate. In this case the epidemiological model turns into a system of ODEs and the condition that host one has smaller case fatality proportion takes the form
where
is the disease-induced death rate of host i,
is the recovery rate and
is the natural death rate. We surmise that in the general case the characteristic equation for host one may exhibit Hopf bifurcation and have complex roots with positive real part, leading to oscillations in the system. However, we were not able to compose an example of that case.
Returning to the immuno-epidemiological model, we consider a continuum of hosts, stratified by their immune response stimulation rate a and the case fatality proportion as a function of a, . We argue that as a function, the case fatality proportion has a minimum because a trade-off exists between the virulence created by the pathogen and the virulence created by the immune response. From the continuum of hosts, the host that has the immune response activation rate
that minimizes
will outcompete the rest. We further find that if the pathogen increases its reproduction rate r, the optimal immune response rate of the host
also increases. This is in agreement with prior results where a different principle for host evolution was used [Citation26]. Finally, we use adaptive dynamics techniques to argue that the minimization of the case fatality proportion is an ESS for the host. Hosts that play nearby strategies cannot invade the host that plays the optimal strategy, that is the host with immune response activation rate
.
We assume logistic growth for the host as opposed to constant growth typically assumed in the literature. The reason for that is that we are looking at the evolution of the host and the host must have an option to persist and an option to die out if not fit. Constant recruitment does not give these two options; it only gives an option to persist, so evolutionary questions about the host cannot be asked in a constant recruitment model. The disease-induced mortality is assumed as a linear combination of the pathogen load and immune response. This linear combination leads to trade-off between mortality due to pathogen and mortality due to immune response and ultimately to the minimum of the CFP as explained in Section 4. Other forms, such as the one assumed by Gilchrist and Sasaki [Citation17] will also lead to a minimum in the CFP. However, some forms of the mortality rate coupled with the form of may result in monotone CFP. Gilchrist and Sasaki [Citation17] and Pugliese [Citation26] assume the transmission rate
proportional to the pathogen load. Comparison to data [Citation22], however, suggests that the transmission rate is not linearly dependent on the pathogen load. Other forms of the transmission rate may presumably also depend on the antibody levels, although the specific form of dependence will have to be determined from experiments and data fitting. At this point we do not expect that adding dependence on the antibody levels will change significantly the results. However, transmission rates that incorporate different susceptibilities of the hosts will lead to a criterion analogous to the one derived by Bowers [Citation9].
One interesting question that remains is whether the principle will remain if we increase the biological complexity of the model. For instance, is the principle valid in vector-host infectious disease-models or environmental transmission models? Our preliminary computations with an immuno-epidemiological vector-host model of arboviral diseases suggest that the principle remains valid [Citation18].
In conclusion, the principle for host evolution that we derive, namely minimization of the case fatality proportion, extends prior results in [Citation17,Citation26] where maximization of host lifespan in the infected class was considered. Our principle is more general in the sense that it applies both to chronic diseases and to diseases with recovery.
Acknowledgments
This project may have benefitted from discussions with Hayriye Gulbudak and Vincent Cannataro. Authors also thank two referees for their helpful suggestions.
Disclosure statement
No potential conflict of interest was reported by the authors.
ORCID
Necibe Tuncer http://orcid.org/0000-0002-6388-2499
Additional information
Funding
References
- S. Alizon and S. Lion, Within-host parasite cooperation and the evolution of virulence, Proc. R. Soc. B 278 (2011), pp. 3738–3747. doi: 10.1098/rspb.2011.0471
- S. Alizon, F. Luciani, and R.R. Regoes, Epidemiological and clinical consequences of within-host evolution, Trends Microbiol. 19(1) (2011), pp. 24–32. doi: 10.1016/j.tim.2010.09.005
- R.M. Anderson and R.M. May, Coevolution of hosts and parasites, Parasitology 85(2) (1982), pp. 411–426. doi: 10.1017/S0031182000055360
- J. -B. André and S. Gandon, Vaccination, within-host dynamics, and virulence evolution, Evolution 60(1) (2006), pp. 13–23. doi: 10.1111/j.0014-3820.2006.tb01077.x
- C.L. Ball, M.A. Gilchrist, and D. Coombs, Modeling within-host evolution of HIV: Mutation, competition and strain replacement, Bull. Math. Biol. 69 (2007), pp. 2361–2385. doi: 10.1007/s11538-007-9223-z
- L.B. Barreiro and L. Quintana-Murci, From evolutionary genetics to human immunology: How selection shapes host defence genes, Nat. Rev. Genet. 11 (2010), pp. 17–30. doi: 10.1038/nrg2698
- S. Bonhoeffer, R.E. Lenski, and D. Ebert, The curse of the pharaoh: The evolution of virulence in pathogens with long living propagules, Proc. R. Soc. London B, Biol. Sci. 263 (1996), pp. 715–721. doi: 10.1098/rspb.1996.0107
- M. Boots and Y. Haraguchi, The evolution of costly resistance in host-parasite systems, Am. Nat. 153 (1999), pp. 359–370. doi: 10.1086/303181
- R.G. Bowers, The basic depression ratio of the host: the evolution of host resistance to microparasites, Proc. R. Soc. Lond. B 268 (2001), pp. 243–250. doi: 10.1098/rspb.2000.1360
- A. Brännström, J. Johansson, and N. von Festenberg, The hitchhiker's guide to adaptive dynamics, Games 4(3) (2013), pp. 304–328. doi: 10.3390/g4030304
- H.J. Bremermann and J. Pickering, A game-theoretical model of parasite virulence, J. Theoret. Biol. 100 (1983), pp. 411–426. doi: 10.1016/0022-5193(83)90438-1
- H.J. Bremermann and H.R. Thieme, A competitive exclusion principle for pathogen virulence, J. Math. Biol. 27 (1989), pp. 179–190. doi: 10.1007/BF00276102
- D. Claessen and A.M. de Roos, Evolution of virulence in a host–pathogen system with local pathogen transmission, Oikos 74 (1995), pp. 401–413. doi: 10.2307/3545985
- D. Coombs, M.A. Gilchrist, and C.L. Ball, Evaluating the importance of within- and between-host selection pressures on the evolution of chronic pathogens, Theor. Pop. Biol 72 (2007), pp. 576–591. doi: 10.1016/j.tpb.2007.08.005
- S. Debroy and M. Martcheva, Immuno-epidemiology and HIV/AIDS: A modeling prospective, in Mathematical Biology Research Trends, Lachlan B. Wilson, ed., Nova Publishers, New York, 2008, pp. 175–192.
- A. Giafis and R.G. Bowers, The adaptive dynamics of the evolution of host resistance to indirectly transmitted microparasites, Math. Biosci. 210 (2007), pp. 668–679. doi: 10.1016/j.mbs.2007.07.006
- M.A. Gilchrist and A. Sasaki, Modeling host-parasite coevolution: A nested approach based on mechanistic models, J. Theor. Biol. 218(3) (2002), pp. 289–308. doi: 10.1006/jtbi.2002.3076
- H. Gulbudak, V. Cannataro, N. Tuncer, and M. Martcheva, Modeling host-parasite coevolution in a coupled arbovirus within-host and between-host system, (in preparation).
- P. Lemey, A. Rambaut, and O. Pybus, HIV evolutionary dynamics within and among hosts, AIDS Rev. 8 (2006), pp. 125–140.
- M. Lipsitch and M.A. Nowak, The evolution of virulence in sexually transmitted HIV/AIDS, J. Theoret. Biol. 174 (1995), pp. 427–440. doi: 10.1006/jtbi.1995.0109
- M. Lipsitch, S. Siller, and M.A. Nowak, The evolution of virulence in pathogens with vertical and horizontal transmission, Evolution 50 (1996), pp. 1729–1741. doi: 10.2307/2410731
- M. Martcheva, An Introduction to Mathematical Epidemiology, Springer, New York, 2015.
- M. Martcheva and X. -Z. Li, Linking immunological and epidemiological dynamics of HIV: The case of super-infection, J. Biol. Dyn. 7(1) (2013), pp. 161–182. doi: 10.1080/17513758.2013.820358
- M.R. Miller, A. White, and M. Boots, The evolution of host resistance: Tolerance and control as distinct strategies, J. Theor. Biol 236 (2005), pp. 198–207. doi: 10.1016/j.jtbi.2005.03.005
- M.A. Nowak and R.M. May, Superinfection and the evolution of parasite virulence, Proc. R. Soc. London, B 255 (1994), pp. 81–89. doi: 10.1098/rspb.1994.0012
- A. Pugliese, The role of host population heterogeneity in the evolution of virulence, J. Biol. Dyn. 5 (2011), pp. 104–119. doi: 10.1080/17513758.2010.519404
- H.R. Thieme and C. Castillo-Chavez, How may infection-age-dependent infectivity affect the dynamics of HIV/AIDS?, SIAM J. Appl. Math. 53 (1993), pp. 1447–1479. doi: 10.1137/0153068
- C.A. Toft and A.J. Karter, Parasite-host coevolution, Trends Ecol. Evol. 5(10) (1990), pp. 326–329. doi: 10.1016/0169-5347(90)90179-H