
THE DISCOUNTED REPRODUCTIVE NUMBER FOR EPIDEMIOLOGY
Timothy C Reluga
Jan Medlock
Alison Galvani
Abstract
The basic reproductive number,, and the effective reproductive number,
, are commonly used in mathematical epidemiology as summary statistics for the size and controllability of epidemics. However, these commonly used reproductive numbers can be misleading when applied to predict pathogen evolution because they do not incorporate the impact of the timing of events in the life-history cycle of the pathogen. To study evolution problems where the host population size is changing, measures like the ultimate proliferation rate must be used. A third measure of reproductive success, which combines properties of both the basic reproductive number and the ultimate proliferation rate, is the discounted reproductive number
. The discounted reproductive number is a measure of reproductive success that is an individual’s expected lifetime offspring production discounted by the background population growth rate. Here, we draw attention to the discounted reproductive number by providing an explicit definition and a systematic application framework. We describe how the discounted reproductive number overcomes the limitations of both the standard reproductive numbers and proliferation rates, and show that
is closely connected to Fisher’s reproductive values for different life-history stages
Key words and phrases: reproductive number, ultimate proliferation rate, game theory
1. Introduction
The theory of biological evolution provides a dynamic method for predicting past and future population change. If we know a population’s current composition and we can identify the fittest individuals, we can predict that population’s future composition. However, when converting theory to prediction we stumble at one particular point: how does nature determine fitness? Fitness is a flexible concept that can be defined in different ways depending on the evolutionary problem being addressed (Mayr, 1997). The evolutionary ecology concept of fitness is conventionally defined as the “reproductive success” of a phenotype or adaptation. This convention bypasses the complexities of population genetics (that inherited genotypes differ from expressed phenotypes according to Mendel’s laws) but is a reasonable approximation under weak selection (Lande, 1982;Charlesworth, 1990).
In infectious disease research, the most commonly invoked measure of reproductive success is the basic reproductive number (Dublin and Lotka, 1925). In a demographic sense,
is the total number of offspring a typical individual expects to have over the course of a lifetime in a naive population. For infectious diseases, this means the number of new infections caused by a typical case in a completely susceptible population. In several archetypal epidemiology models, strains of disease with higher
values out-compete strains with lower
values, leading to a rule-of-thumb that evolution maximizes
. This phenomenon has been exploited to derive a number of useful results concerning pathogen evolution in the context of vaccination (Gandon et al., 2001,2003). In populations that are not completely susceptible,
is replaced by an effective reproductive number
representing the number of new cases created, conditional on the population’s state. ButBremermann and Thieme (1989) have shown that even in some cases where populations are not completely susceptible,
can be used to predict the asymptotic outcomes of some nonlinear competition models.
maximization does not always correctly predict evolutionary changes in strain frequencies, however. If, for example, at the beginning of an epidemic, there are two strains of a virus: one strain that kills its host in 2 days while causing 2 new cases at the end of those 2 days, and a second strain that kills its host in 6 days while causing 3 new cases at the end of those 6 days, the second strain has a larger
(3 versus 2). However, the relative frequency of the second strain will decrease over time. In 6 days, the first strain will have been transmitted though 3 generations, for a total of 8 descendant infections compared to the 3 descendant infections of the second strain. If both strains start with a single case, then after a month, there will clearly be more cases of strain one than strain two although strain two has the larger
value.
The above argument applies for both naive populations with and non-naive population states where we use the effective reproductive number
. Both
and
ignore the timing of reproduction events, treating all reproduction events as having the same value. If the population size is not changing, all reproductive events do indeed have the same value, and reproductive success can be predicted by
or
. But if the population size is changing, early reproductive events have more value than late reproduction events.
In situations where the population size is growing, reproductive success is often measured in terms of the ultimate proliferation rater of the strains. Biologically,r is the asymptotic instantaneous rate of change in abundance of a population of identical individuals after the population structure has settled down to a stable stage distribution (Keyfitz, 1968). Although no population can grow indefinitely, and the allowance for indefinite growth superficially violates one of the fundamental premises of Darwin and Wallace’s theory of evolution,r is a measure of reproductive success that receives regular use in both theory and practice (Bulmer, 1994;McGraw and Caswell, 1996;Caswell, 2001). Comparative analysis of the ultimate proliferation ratesr of different strains can correctly predict evolutionary changes in strain frequency not just when population sizes are changing, but also when the total population size is stationary (r = 0).
However, the dichotomy between the ultimate proliferation rate and the basic reproductive value is unsettling. A better mathematical theory of fitness would not require this ad-hoc distinction between populations that are changing size and populations that are not changing size. The ultimate proliferation rate is often more difficult to calculate than a reproductive number when comparing novel mutant strains with a wild-type population, especially when considering a subpopulation having a different life history than the population as a whole. Another issue is that the ultimate proliferation rate has the connotation of a population-scale concept, describing how the size of a specific subpopulation changes in time as the population structure converges to a stationary distribution. By contrast, reproductive numbers have a connotation associated with individuals. As each individual case of an infection has a basic reproductive number, so each mutation in an evolving population has the potential to create a new strain with its own reproductive number.
Instead of having to deal with both the ultimate proliferation rate and the basic reproductive number, there is an equivalent alternative approach. Here we describe a third measure of reproductive success equally valid in contexts of both growing and static background populations: the discounted reproductive number. The discounted reproductive number
is a measure of reproductive success that is an individual’s expected lifetime offspring production discounted by the background population growth rate. We provide an explicit definition and systematic application framework for the discounted reproductive number. Although other authors have used the discounted reproductive number for life-history optimization (Goodman, 1982;McNamara, 1991;McNamara et al., 2001), its merits remain under-recognized. Furthermore, the fundamentally comparative and game-theoretic nature of the discounted reproductive number has not been made explicit. In this work, we illustrate the versatility of the discounted reproductive number as a fitness measure that subsumes the two most common measures, the basic reproductive number and the ultimate proliferation rate. We also show that
is closely connected to Fisher’s reproductive values for different life-history stages. Our results show how the standard concept of reproductive numbers extends naturally to solve a wide array of evolution problems.
2. Mathematical methods and results
To define measures of reproductive success within a quantitative framework, we consider the family of continuous-time, stage-structured population models described by a (possibly nonlinear) matrix equation of the form
(1) |
wheren is a vector of population abundances at each life stage,F is a non-negative matrix of reproduction rates, andT is the negative of an M-matrix1 describing the transition rates of a Markov process, andF +T is irreducible and primitive. Formally, there are many different splittings of the system into fecundity and transition components such thatF +T =F̃ +T̃ withF ≠F̃ andT ≠T̃. Biologists may argue that life-cycle structure dictates a canonical splitting. Our results will hold for all splittings provided the non-degeneracy conditionλ0(F +T)> λ0(T) holds, where the notationλ0(·) represents the matrix’s largest real eigenvalue.
The ultimate proliferation rate is the largest real eigenvalue of the growth operator:
(2) |
which is the eigenvalue with largest real part. The general methods for calculation of and
are closely related to the theory of Markov decision processes developed byHoward (1960) and to branching process theory (Ulam, 1990). The basic reproductive number can be thought of as the expected number of offspring produced in each life-history stage, weighted by the probability of surviving and entering that stage, and summed over all stages. The basic reproduction number of an epidemic is the expected number of new cases in a naive population.
is calculated (Cushing and Yicang, 1994;Caswell, 2001;Diekmann et al., 1990;van den Driessche and Watmough, 2002) as
(3) |
In situations where the expected number of new cases is conditional on a population state other than the naive state, the calculation is the same and the result is called an effective reproductive number. Note that the difference between
and
is simply whether the matricesF andT, as functions of the population state, are evaluated at a naive or non-naive state. In what follows, we will refer only to
, although the statements also hold for
.
A third measure of reproductive success that combines both the basic reproductive number and the ultimate proliferation rate is the discounted reproductive number. The discounted reproductive number represents the number of off-spring an individual expects to have over the course of a lifetime, discounted for expected changes in population size. The conceptual origins of the discounted reproductive number begin with ther–K theory of selection introduced byMargalef (1959) andMacArthur and Wilson (1967). Ther–K selection theory was originally envisioned as a continuum betweenr selected species in resource-rich environments that evolve to maximize productivity andK selected species in resource-scarce environments that evolve to maximize efficiency. The nomenclature derives from the simple logistic-growth model
(4) |
wherer represents the per-capita growth rate andK denotes the carrying capacity. In application, however,r–K selection has most often been treated as a polar dichotomy because there is no natural continuum between the carrying capacityK, measured in the same units as the population’s size, and the growth rater, measured in units of inverse time. Although it is a convenient caricature,r–K selection theory has been largely abandoned.
The discounted reproductive number supplies the missing spectrum tor–K selection theory, withK selection corresponding to static populations without discounting, andr selection corresponding to any scenario with positive population growth. The concept of discounting has a long history in economics tracing back to Fibonacci’s 13th-century work (Goetzmann, 2005). In the later part of the 20th century, some researchers have observed that the discount rate also plays an important role in evolutionary population biology because it provides a means to account for costs associated with delays in reproduction (Goodman, 1982). has appeared repeatedly throughout the historical development of evolutionary theory, although not explicitly by this name.Taylor et al. (1974) showed that the proliferation rate and
are equivalent fitness measures in the special case of a discrete-time, age-structured model.Schaffer (1974) also implicitly employed the discounted reproductive number in his study of life-history evolution.Li and Schneider (2002) later extended this by studying models with general stage structure in discrete time.Goodman (1982) seems to be the first author to treat the discount rate as a free parameter, although he takes pains to downplay any originality in his paper. Goodman further suggests a method for optimizing
in an age-structured model. The method implicitly requires that an equilibrium strategy have a form of convergent stability (Eshel, 1983) to guarantee that it will be found. McNamara refers specifically to
as the “discounted maximum future expected reward” (McNamara, 1991) or the “discounted lifetime reproduction success” (McNamara et al., 2001).
plays an important role in studies of evolution based on dynamic programming approaches to fitness (McNamara, 1991) and alternative versions of adaptive topographies (McNamara, 1993;Lande, 1982).
The discounted reproductive number is defined in a manner similar to the basic reproductive number, except that the number of offspring from each stage is also weighted to reduce the value of offspring produced later in time, relative to offspring produced earlier. Thus, the discounted reproductive number is defined as
(5) |
whereδ is the discount rate. In general, we takeδ > λ0(T) so that an inverse exists. Since all known species are mortal,λ0(T)< 0.
Just like,
can be defined for other common life-history models in addition to the continuous-time matrix models presented here. For a McKendrick–von Foerster model,
(6) |
wherea is the age of a cohort, the discounted reproductive number is given by
(7a) |
with
(7b) |
For the discrete-time matrix model
(8) |
the discounted reproductive number is given by
(9) |
whereθ is the discrete-time discounting ratio (Li and Schneider, 2002). The equation = 1 is a special case which leads to the well-known Euler–Lotka eigen-value problem (Kot, 2001) and Fisher’s reproductive value (Fisher, 1930) for age-structured populations.
3. Properties of
Inspectingequations (3) and(5), we see that when there is no discounting (δ = 0), =
. Further, positive discount rates always result in discounted reproductive numbers smaller than the basic reproductive number.
Proposition 1
The discounted reproductive number is a decreasing function of the discount rateδ and is bounded below by 0.
Proof
Define the discounted generation matrix
(10) |
so that =λ0(G). By the properties of M-matrices,G exists for allδ > λ0(T). The discounted generation matrix can be rewritten as
(11) |
Since multiplication by a scalar commutes over all matrices, we can decompose the matrix exponential to the form
(12) |
Because −T is an M-matrix, the matrix exponential and hence the whole integrand must be non-negative (Berman and Plemmons, 1979, p. 146). Perron–Frobenius theory tells us that the dominant eigenvalue of a non-negative matrix is non-negative (Horn and Johnson, 1985), so =λ0(G) ≥ 0.
By inspection, the integrand is decreasing inδ. BecauseF is also non-negative, we see thatG is a component-wise strictly decreasing function ofδ in the sense that either ∂Gij/∂δ < 0 orGij = 0. Takeh such that for sufficiently smallε,
(13) |
If for two matricesA andB whereB ≥A ≥ 0 component-wise, thenλ0(B) ≥λ0(A) ≥ 0 (Horn and Johnson, 1985, p. 491), so
(14) |
It follows that eitherλ0(G(δ)) = 0, orλ0(G(δ))> λ0(G(δ +ε)) ≥ 0 in the limit asε approaches 0 from above.
As =λ0(G), we conclude that
is decreasing in the discount rateδ, and strictly decreasing so long as
is positive.
Note that this proof holds for anyT andF satisfying our assumptions, so there is a continuous cone of discounted reproductive numbers for any alternative splittingF̃ +T̃ =F +T, but all of these
’s have the same relationship to the ultimate proliferation rater: the discounted reproductive number is 1 if reproduction is discounted at the ultimate proliferation rate.
Proposition 2
= 1 if and only ifδ =r.
Proof
We need to show that the dominant eigenvalue ofF(rI −T)−1 is 1. Sincer is the dominant eigenvalue of the primitive matrixF +T, the left eigenvectorv corresponding tor is unique (up to scaling) and can be taken to be positive (Horn and Johnson, 1985). Moreover,v is the only positive left eigenvector for any eigenvalue (Chow, 1983, p. 914). Then
(15) |
wherevT denotes the transpose ofv. Since −T is an M-matrix andr =λ0(F+T)> λ0(T), the matrixrI −T also satisfies the definition of an M-matrix. Because the inverse of an M-matrix is non-negative (Horn and Johnson, 1991, p. 117), (rI −T)−1 exists and is non-negative. If we multiply by this inverse,
(16) |
We now see thatv is a positive left eigenvector ofF(rI −T)−1 corresponding to the eigenvalue of 1.
Perron–Frobenius theory tells us that for a non-negative matrixA, if an eigenvaluee has a modulus smaller than the dominant modulus, then the eigenvector fore is not strictly positive. Contrapositively, ifA is non-negative andv is a positive eigenvector, then the corresponding eigenvalue has modulus equal to the dominant eigenvalue (Chow, 1983, p. 914). Sincev is a positive left eigenvector of the nonnegative matrixF(rI −T)−1, we must haveλ0(F(rI −T)−1) = 1. Thus, ifδ =r, = 1.
Conversely, letδ ≠r. By Proposition 1, is strictly decreasing inδ when
> 0, so
≠ 1.
4. Comparison of fitness measures
To demonstrate that the discounted reproductive number is equivalent to the basic reproductive number and the ultimate proliferation rate as a fitness measure, consider a very large population where a small number of individuals possess a mutant phenotypeπ that differs from the population’s resident phenotypeπ̄. The life-history and fecundity matrices of individuals depend both on the mutant phenotype and on the resident phenotype, which we will denote using subscripts. The mutant phenotype has basic reproductive number
(17) |
ultimate proliferation rate
(18) |
and discounted reproductive number
(19) |
where
(20) |
We have chosen the discount rateδ(π̄) to be the ultimate proliferation rate of individuals with the resident phenotype because this is the population’s baseline growth rate against which the mutant phenotype is competing. Although all phenotypes with(π,π̄)> 1 will increase in number, only phenotypes with an ultimate proliferation rate larger than that of the resident population will increase in frequency.
A mutant phenotype can only invade a population if it has a greater fitness in the current environment than the resident phenotype. If we define the relative fitness of a phenotype asW(π,π̄), with the convention thatW(π̄,π̄) = 1, then a mutation can invade whenW(π,π̄)> 1. A resident phenotype is a strict Nash equilibrium of the population game if all possible mutations have a lower fitness than the resident phenotype (W(π,π̄)< 1 for allπ ≠π̄). This and other evolutionary properties of a fitness measure can be assessed using the pairwise invasion test map
(21) |
As examples, we re-derive and generalize results fromMylius and Diekmann (1995) andBulmer (1994, p. 77) describing the circumstances when each fitness measure is equivalent.
Comparison 1
When a population is at a stable size (i.e., no growth,r = 0), the pairwise invasion test of relative fitness measured in terms of the basic reproductive number is the same as the pairwise invasion test of relative fitness measured in terms of the discounted reproductive number.
Proof
We need to prove
(22) |
Because there is no asymptotic growth of the population (r(π̄,π̄) = 0), there is no discounting (δ(π̄) = 0) and the discounted reproductive numbers are equal to the basic reproductive numbers. The ratio of the mutant discounted reproductive number to the resident discounted reproductive number is then the same as the ratio of the basic reproductive numbers:
(23) |
Thus, the population games are identical, and the pairwise invasion tests must also be identical.
does not correctly measure fitness, however, when a population’s size is growing or decaying because it gives all reproduction events equal weight. In a growing population, reproduction events next year are worth less than reproduction events this year. Both the ultimate proliferation rate and the discounted reproductive number are suitable alternatives in these cases.
Comparison 2
The pairwise invasion test of relative fitness measured in terms of the ultimate proliferation rate is the same as the pairwise invasion test of relative fitness measured in terms of the discounted reproductive number.
Proof
In order for 1 to remain the critical value for invasion, as in the definition (21) ofS(W), we will use the exponential of the ultimate proliferation rate, er(π,π̄), instead of the ultimate proliferation rate. Consequently, the relative fitness ofπ with respect toπ̄ is
(24) |
Thus, we need to show that
(25) |
The tests are equal if they are equal at each point (π,π̄). Ifr(π,π̄) −r(π̄,π̄) = 0, then the discounted reproductive number has a discount rate equal to the proliferation rate and by Proposition 2,(π,π̄) = 1. Now by Proposition 1, ifr(π,π̄)−r(π̄,π̄)< 0, then the discount rate is faster than the proliferation rate and
(π,π̄)< 1. Ifr(π,π̄) −r(π̄,π̄)> 0, then the discount rate is slower than the proliferation rate and
(π,π̄)> 1. The tests are therefore equal.
Comparisons 1 and 2 show that the discounted reproductive number subsumes the two most common measures of fitness,r and. To show how
can be combined with population-level and individual-level models to study evolutionary strategies, letn(t) be the macroscopic state of the resident population,e(t) be the state of the environment, andp(t) be the probability density of the life-history state space over time for an individual with the mutant phenotypeπ. We refer top(t) as the microscopic state, because it is the life-history state of a single individual. We refer ton(t) as the population’s macroscopic state because it is a combination of all other individuals’ states. Mathematically, the macroscopic and microscopic dynamics are driven by separate but related processes. From the microscopic perspective of the individual, life is stochastic, but from the macroscopic perspective of a large population, the dynamics are treated as deterministic. If we use a Markov process to describe the life of an individual, and a system of ordinary differential equations to describe changes in a population, the relative fitness of an individual with phenotypeπ is given by the discounted reproductive number
(26a) |
where
(26b) |
governs the dynamics of the microscopic state, and
(26c) |
(26d) |
govern the dynamics of the macroscopic state and the environment. Here,T is the transition-rate matrix of the life-history Markov process,δ =r(π̄,π̄) is the population’s proliferation rate,F is the fecundity matrix, and the vectorv is the reproductive value for each state. The use ofv was introduced byFisher (1930), who suggested it be interpreted as the relative expected future lifetime contribution to reproductive output for individuals in each state whenπ =π̄. The population’s proliferation rateδ discounts future reproduction relative to current reproduction.
Consider a system that is stationary in the sense that the resident population is growing at rateδ(π̄) with fixed proportions of life-history states so thatn(t) = eδ(π̄)tn* and the resource is at steady state, i.e.e(t) =e*. The dynamic equations can be reduced by substitutingn(t) = eδ(π̄)tn* ande(t) =e* into system (26) and requiring that the equation for the microscopic state,p, be at steady state, resulting in the system of equations:
(27a) |
(27b) |
(27c) |
Under these equilibrium conditions, the probability of being in each state at a given time is
(28) |
It follows that2
(30) |
We want to choose the reproductive valuesv such that the relative fitnessW(π,π̄) correctly predicts changes in strategy frequency of long times; mathematically,
(31) |
If we choosev to be the left eigenvector associated with the dominant eigenvalue ofF (δI −T)−1, normalized such thatvTp(0) = 1, then
(32) |
Applying Comparison 2, we find that for this choice ofv,
(33) |
Thus, the relative fitnessW(π,π̄) correctly predicts which strategiesπ can replace the resident strategyπ̄ provided that the reproductive valuev is defined in terms of a left eigenvector. The reproductive valuev is a generalization of Fisher’s reproductive value because it depends on both the invading phenotypeπ and the resident phenotypeπ̄, whereas for Fisher’s reproductive value it is traditionally assumed thatπ =π̄. The definition of the reproductive values in terms of an eigenvector ofF (δI −T)−1 may seem to preclude the practical use ofequation (26a). However, oftenF has only 1 nonzero row, corresponding to birth into only one life-history stage, andvTF is a scalar multiple of that row. In this case, the conditionvTp(0) = 1 completely characterizesvT.
Given a fitness function, we can then compare strategies and determine which strategies are Nash equilibria. A phenotypeπ̄ is a Nash equilibrium strategy if(π*,π*)>
(π*,π*) for allπ ≠π. When
is a smooth function, the Nash condition is equivalent to and.
5. Applications of
The discounted reproductive number is easy to apply, although closed-form solutions are uncommon. As a first example, let us calculate for a simplification of the classic problem of maturation-time evolution. Consider a population with two life-history states, juvenile and adult, that grows linearly according to
(34) |
(35) |
whereπ̄ is the resident phenotype strategy. An individual with phenotypeπ has life-history transition-rate matrix
(36) |
and fecundity matrix
(37) |
Applyingequation (27), we find the discounted reproductive number
(38) |
where the discount rate
(39) |
Let us extend this model to pathogen evolution in epidemiology. Consider incubation-period evolution in the susceptible-exposed-infectious-recovered (SEIR) model
(40a) |
(40b) |
(40c) |
(40d) |
HereS is the number of people who are susceptible,E is the number exposed,I is the number infectious, andR is the number recovered. The birth rate is given byη,m is the background death rate,μE is the infection-induced death rate during the exposed stage,μI is the infection-induced death rate during the infectious stage, andγ is the recovery rate. The infection rate,β, and the incubation rate,f, are both functions of pathogen phenotype.
The dynamics of System (40) are well understood. If
(41) |
then the disease-free stationary solution is globally attracting. Otherwise, there is a unique, globally attracting stationary solution with known Lyapunov function (Korobeinikov and Maini, 2005).
In infectious-disease modeling, the pathogens are the organisms of interest, and they only exist in infected human hosts. Therefore, we only need the exposed (E) and infectious (I) compartments of the host population to model the pathogen population. The life-history transition-rate matrix for an infection with typeπ is
(42) |
and the fecundity matrix for the creation of new infections is
(43) |
We applyequation (5) and find
(44) |
with
(45) |
In this calculation,S is treated as a constant. There are two natural cases ofequation (40) where this is reasonable. First, early in an epidemic, the number of infections is small whileS is large and relatively constant. In this case,S is fixed and independent of the resident phenotypeπ̄. Second, late in an epidemic, the dynamics approach the endemic steady-state solution, whereS is again constant. In this case, the steady-state depends on the resident phenotypeπ̄;
(46) |
As an example, letf(π) =π and (Fig. 1). Example contour plots of(π, π̄) early and late in an epidemic are shown inFigure 2. The condition above are used to show that the points where the 1-contours cross inFigure 2 are Nash equilibria.
Figure 1.
Transmission rate,β vs. incubation rate,f used in our example forequation (40) andFigure 2. The functional forms aref(π) =π and withb0 =b1 = 1.
Figure 2.
Contour plots of(π, π̄) early (top) and late (bottom) in the course of an epidemic described byequation (40). Nash equilibria correspond to the phenotypes where the 1-contours intersect. At the start of an epidemic, the number of susceptibles is independent of the strategies and the number of infections is growing exponentially. However, near the end of an epidemic, the number of susceptible individuals at dynamic equilibrium depends on the population’s resident phenotypeπ̄, although there is no net change in the number of susceptible individuals per unit time. Because of these differences, the Nash equilibriumπ* for the rate of leaving the exposed class is more than twice as large early compared to late in the epidemic. Parameter valuesμE =μI = 0,γ = 10,m = 1/60,η/m = 100. The early plot is for the naive population,S = 100. For the late plot,S is given byequation (46).
The value of the discount rate differs at different points in an epidemic, as does the number of susceptiblesS. Early in an epidemic almost all individuals are susceptible and the number of infections grows exponentially. Thus, the discount rate is large. Over time, however, the excess of susceptibles is exhausted and the population-dynamics approach an endemic stationary solution with the discount rate also approaching 0. The evolutionary pressures on the pathogen change as the discount rate varies between these extremes. The discount rate is faster early in the epidemic when the number of susceptibles is largest and the number of infections is growing quickly, favoring short incubation periods (largeπ). Late in the epidemic, discounting is slow, and longer incubation periods can pay off with greater transmission.Figure 2 shows how these differences play out in moving the Nash equilibrium from short incubation early in the epidemic to longer incubation periods late in the epidemic.
6. Discussion
We have shown that the discounted reproductive number unifies several disparate measures of fitness within a single framework while providing a mechanistic basis for the derivation of new results. The discounted reproductive number allows us to systematically account for the timing of transmission events to study changes in frequency in a way not possible with the reproductive number
. Although we illustrated the use of
with ordinary differential equation models, we also showed how
is formulated for difference equation and partial differential equation models. We have applied
to model changes in host behavior when virulence changes with the age of the host (Reluga et al., 2007) and here to describe changes in pathogen selection over the course of an epidemic. In addition,
can be used in studies of the evolution of clutch size, senescence, growth, and other evolutionary trade-offs, including co-evolution, like that between host and pathogen. We also have shown that
(π,π̄) is closely connected to Fisher’s reproductive values for different life-history strategies. Therefore, the discounted reproductive number is broadly applicable to other studies where it is reasonable to assume populations are near dynamic equilibrium or growing exponentially.
For many models, the reproductive number is easier to derive and takes a simpler form than the ultimate proliferation rater. When discount rateδ is treated as a free parameter, which may suffice for some kinds of analysis,
is as easy to derive and as simple in form as
. However, whenδ is taken to be the ultimate proliferation rate of the resident populationr(π̄,π̄), the determination of
(π,π̄) involves the complexity of calculating bothr and
.
The proofs of Proposition 1 and Proposition 2 rely on our assumption that the demography matrixF +T is irreducible and primitive. It is not immediately clear what form the results would take if either of these conditions were relaxed. Of particular interest is the fact that the splitting ofF +T needed to define, as well as
, is not unique. A similar issue is present and more onerous when matrix models are generalized to branching-process models of demography. This suggests two points. First, the mathematics is not capturing some parts of the biology, since fecundity can not be defined arbitrarily in practice. Second, while each splitting defines a different
, there seem to be no significant mathematical differences between these splittings in practice. These points suggest a need for a practical convention to aid in comparative analyses. Some readers may note a close connection between the concept of the discounted reproductive number as a function of the discount rateδ and the concept of a Laplace transform. This connection may provide a convenient avenue for mathematical applications of
in new research.
Evolutionary population game theory based on the discounted reproductive number has a number of limitations. While convenient, population game theory based on is an asymptotic theory with well-known shortcomings (Abrams et al., 1993). The theory has primarily been applied in settings where population dynamics are stationary or growing exponentially.Equation (5) assumes the matricesF andT are stationary. If transition rates or fecundities are not stationary, the integrals inequation (11) andequation (26a) do not have general closed-form solutions. Similar issues arise in attempts to apply
andr in situations where dynamics oscillate. Further generalizations are needed to accommodate dynamic programming formulations and Lyapunov-exponent descriptions of fitness in settings with periodic and chaotic population dynamics (McNamara et al., 2001 ;Metz et al., 1992 ). We have also sidestepped the challenges associated with frequency-dependent selection by restricting our attention to population games. In general, alternative mathematical approaches should be used if the evolutionary dynamics of a biological system are not easily summarized in terms of asymptotic properties.
In summary, we have provided an explicit definition and systematic application framework for extending the reproductive number to the discounted reproductive number. We demonstrated that the discounted reproductive number is a versatile fitness measure that subsumes the two most common measures of fitness, the basic reproductive number and the ultimate proliferation rate.
Figure 3.
Nash equilibrium strategy (π*) at the start of an epidemic, depending on the initial population size (N =η/m). Parameter values arem = 1/60,μE =μI = 0,γ = 1/3, andb0 andb1 as inFigure 1.
Acknowledgments
The authors thank the Notsew Orm Sands Foundation for financial support. Portions of this work were performed under the auspices of the U.S. Department of Energy under contract DE-AC52-06NA25396. This work was supported in part by NIH grants AI28433 and RR06555 (ASP) and the Human Frontiers Science Program grant RPG0010/2004. JM was supported by NIH grant 2 T32 MH020031-07 for this work.
This manuscript is dedicated to Fred Brauer and Karl Had-eler, for their valuable contributions to the field of mathematical epidemiology. T. Reluga thanks Gardner Brown for initiating discussions.
Terminology
Here we briefly outline the terminology we use to describe a variety of matrices that arise in population modeling.
A matrixA = [aij] isnon-negative if all its entries are non-negative,aij ≥ 0. Likewise,A ispositive if its entries are all positive,aij > 0.
AZ-matrix is a square matrix whose off-diagonal entries are all non-positive,aij ≤ 0 fori ≠j. A Z-matrix can be written asZ =αI −P, whereP is a non-negative matrix. An important property of Z-matrices is that the matrix exponential e−Zt is non-negative for allt ≥ 0 (Berman and Plemmons, 1979, p. 146).
AnM-matrix is a Z-matrix that can be decomposed asM =αI −P for some non-negative matrixP and someα > λ0(P), whereλ0 denotes the largest real eigenvalue of the matrix. The real part of all the eigenvalues of an M-matrix are positive. Moreover, the inverse of an M-matrix is non-negative (Horn and Johnson, 1991, p. 117). SeeHorn and Johnson ( 1991) andBerman and Plemmons ( 1979) for expositions on M-matrix theory.
A matrixA beingirreducible means, in population terms, that each life stage can ultimately produce every other life stage, either by transition between life stages or by production of offspring. (See, for example,Caswell ( 2001, p. 81).) Note that this definition implies that the matrix must be square. An example of a reducible matrix would be the projection matrix from a model that includes post-reproductive adults: the post-reproductive life stage cannot produce newborns or any other life stage.
A non-negative matrixP isprimitive if there exists an integerz such thatPz is positive (Horn and Johnson, 1985). Note that primitive implies irreducible. Otherwise, we define the matrixA to beprimitive if there existsα > 0 such thatA =P −αI, withP non-negative and primitive. Note thatA primitive implies that it is irreducible, as with non-negative matrices, and that it is the negative of a Z-matrix,A = −Z. For a primitive matrixA =P −αI, by Perron–Frobenius theory, the eigenvalue of maximum modulus ofP is real and positive (Horn and Johnson, 1985). Therefore, the largest real eigenvalue,λ0(P), is the eigenvalue of maximum modulus. This implies thatλ0(A) =λ0(P) −α is real and is the eigenvalue ofA with the largest real part.
Footnotes
See theappendix for a definition of M-matrix and other matrix terminology.
(29) |
Contributor Information
Timothy C. Reluga, Department of Mathematics, Pennsylvania State University, State College, PA 16802
Jan Medlock, Department of Epidemiology and Public Health, Yale University School of Medicine, New Haven, CT 06520.
Alison Galvani, Department of Epidemiology and Public Health, Yale University School of Medicine, New Haven, CT 06520.
References
- Abrams PA, Matsuda H, Harada Y. Evolutionarily unstable fitness maxima and stable fitness minima of continuous traits. Evolutionary Ecology. 1993;7:465–487. [Google Scholar]
- Berman A, Plemmons RJ. Nonnegative Matrices in the Mathematical Sciences. Academic Press; New York, NY: 1979. [Google Scholar]
- Bremermann HJ, Thieme HR. A competitive-exclusion principle for pathogen virulence. Journal of Mathematical Biology. 1989;27:179–190. doi: 10.1007/BF00276102. [DOI] [PubMed] [Google Scholar]
- Bulmer M. Theoretical Evolutionary Ecology. Sinauer; Sunderland, MA: 1994. [Google Scholar]
- Caswell H. Matrix Population Models. 2. Sinauer; Sunderland, MA: 2001. [Google Scholar]
- Charlesworth B. Optimization models, quantitative genetics, and mutation. Evolution. 1990;44:520–538. doi: 10.1111/j.1558-5646.1990.tb05936.x. [DOI] [PubMed] [Google Scholar]
- Chow T-S. Matrices and linear algebra. In: Pearson C, editor. Handbook of Applied Mathematics: Selected Results and Methods. 2. Van Nostrand Reinhold; New York: 1983. pp. 878–927. [Google Scholar]
- Cushing JM, Yicang Z. The net reproduction value and stability in matrix population models. Natural Resources Modeling. 1994;8:297–333. [Google Scholar]
- Diekmann O, Heesterbeek JAP, Metz JAJ. On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. Journal of Mathematical Biology. 1990;28:365–382. doi: 10.1007/BF00178324. [DOI] [PubMed] [Google Scholar]
- Dublin LI, Lotka AJ. On the true rate of natural increase. Journal of the American Statistical Association. 1925;20:305–339. [Google Scholar]
- Eshel I. Evolutionary and continuous stability. Journal of Theoretical Biology. 1983;103:99–111. doi: 10.1006/jtbi.1996.0312. [DOI] [PubMed] [Google Scholar]
- Fisher RA. The Genetical Theory of Natural Selection. Oxford University Press; Oxford, UK: 1930. [Google Scholar]
- Gandon S, Mackinnon M, Nee S, Read A. Imperfect vaccination: some epidemiological and evolutionary consequences. Proceedings of the Royal Society of London Series B-Biological Sciences. 2003;270:1129–1136. doi: 10.1098/rspb.2003.2370. [DOI] [PMC free article] [PubMed] [Google Scholar]
- Gandon S, MacKinnon MJ, Nee S, Read AF. Imperfect vaccines and the evolution of pathogen virulence. Nature. 2001;414:751–756. doi: 10.1038/414751a. [DOI] [PubMed] [Google Scholar]
- Goetzmann WN. Fibonacci and the financial revolution. In: Goetzmann WN, Rouwenhorst KG, editors. The origins of value. The financial innovations that created modern capital markets. Oxford University Press; New York, NY: 2005. pp. 123–144. [Google Scholar]
- Goodman D. Optimal life histories, optimal notation, and the value of reproductive value. American Naturalist. 1982;119:803–823. [Google Scholar]
- Horn RA, Johnson CR. Matrix Analysis. Cambridge University Press; Cambridge, UK: 1985. [Google Scholar]
- Horn RA, Johnson CR. Topics in Matrix Analysis. Cambridge University Press; Cambridge, UK: 1991. [Google Scholar]
- Howard RA. Dynamic Programming and Markov Processes. MIT Press; Cambridge, MA: 1960. [Google Scholar]
- Keyfitz N. Introduction to the Mathematics of Population. Addison-Wesley; Reading, Massachusetts: 1968. [Google Scholar]
- Korobeinikov A, Maini P. Non-linear incidence and stability of infectious disease models. Mathematical Medicine and Biology. 2005;22:113–128. doi: 10.1093/imammb/dqi001. [DOI] [PubMed] [Google Scholar]
- Kot M. Elements of Mathematical Ecology. Cambridge; New York, NY: 2001. [Google Scholar]
- Lande R. A quantitative genetic theory of life history evolution. Ecology. 1982;63:607–615. [Google Scholar]
- Li CK, Schneider H. Applications of Perron-Frobenius theory to population dynamics. Journal of Mathematical Biology. 2002;44:450–462. doi: 10.1007/s002850100132. [DOI] [PubMed] [Google Scholar]
- MacArthur RH, Wilson EO. The Theory of Island Biogeography, Monographs in Population Biology. Princeton University Press; Princeton, NJ: 1967. [Google Scholar]
- Margalef R. Mode of evolution of species in relation to their places in ecological succession. In: Hewer HR, Riley ND, editors. XVth international congress of zoology. 1959. pp. 787–789. [Google Scholar]
- Mayr E. The objects of selection. Proceedings of the National Academy of Sciences of the United States of America. 1997;94:2091–2094. doi: 10.1073/pnas.94.6.2091. [DOI] [PMC free article] [PubMed] [Google Scholar]
- McGraw JB, Caswell H. Estimation of individual fitness from life-history data. The American Naturalist. 1996;147:47–64. [Google Scholar]
- McNamara JM. Optimal life histories — a generalization of the Perron–Frobenius theorem. Theoretical Population Biology. 1991;40:230–245. [Google Scholar]
- McNamara JM. Evolutionary paths in strategy space — an improvement algorithm for life-history strategies. Journal of Theoretical Biology. 1993;161:23–37. doi: 10.1006/jtbi.1993.1037. [DOI] [PubMed] [Google Scholar]
- McNamara JM, Houston AI, Collins EJ. Optimality models in behavioral biology. SIAM Review. 2001;43:413–466. [Google Scholar]
- Metz JAJ, Nisbet RM, Geritz SAH. How should we define fitness for general ecological scenarios. Trends in Ecology and Evolution. 1992;7:198–202. doi: 10.1016/0169-5347(92)90073-K. [DOI] [PubMed] [Google Scholar]
- Mylius SD, Diekmann O. On evolutionarily stable life histories, optimization and the need to be specific about density dependence. Oikos. 1995;74:218–224. [Google Scholar]
- Reluga TC, Medlock J, Poolman E, Galvani AP. Optimal timing of disease transmission in an age-structured population. Bulletin of Mathematical Biology. 2007;69:2711–2722. doi: 10.1007/s11538-007-9238-5. [DOI] [PubMed] [Google Scholar]
- Schaffer WM. Selection for optimal life histories — effects of age structure. Ecology. 1974;55:291–303. [Google Scholar]
- Taylor HM, Gourley RS, Lawrence CE, Kaplan RS. Natural-selection of life-history attributes — analytical approach. Theoretical Population Biology. 1974;5:104–122. doi: 10.1016/0040-5809(74)90053-7. [DOI] [PubMed] [Google Scholar]
- Ulam SM. Los Alamos Series in Basic and Applied Sciences. University of California Press; 1990. Analogies Between Analogies: The mathematical reports of S. M. Ulam and his Los Alamos collaborators. [DOI] [PubMed] [Google Scholar]
- van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences. 2002;180:29–48. doi: 10.1016/s0025-5564(02)00108-6. [DOI] [PubMed] [Google Scholar]