Abstract
The basic reproductive ratio, , is one of the fundamental concepts in mathematical biology. It is a threshold parameter, intended to quantify the spread of disease by estimating the average number of secondary infections in a wholly susceptible population, giving an indication of the invasion strength of an epidemic: if , the disease dies out, whereas if , the disease persists. has been widely used as a measure of disease strength to estimate the effectiveness of control measures and to form the backbone of disease-management policy. However, in almost every aspect that matters, is flawed. Diseases can persist with , while diseases with can die out. We show that the same model of malaria gives many different values of , depending on the method used, with the sole common property that they have a threshold at 1. We also survey estimated values of for a variety of diseases, and examine some of the alternatives that have been proposed. If is to be used, it must be accompanied by caveats about the method of calculation, underlying model assumptions and evidence that it is actually a threshold. Otherwise, the concept is meaningless.
1. Introduction
The basic reproductive ratioβalso known as the basic reproductive number, the basic reproduction number, the control reproduction number, or βis one of the foremost concepts in epidemiology [1β3]. is the most widely used epidemiological measurement of the transmission potential in a given population [4]. It is a measure of initial disease spread, such that if , then the disease can invade an otherwise susceptible population and hence persist, whereas if , the disease cannot successfully invade and will die out. The concept is defined as the number of secondary infections produced by a single infectious individual in an otherwise susceptible population [5].
Despite its place at the forefront of mathematical epidemiology, the concept of is deeply flawed. Defining proves to be significantly more difficult than it appears. Few epidemics are ever observed at the moment an infected individual enters a susceptible population, so calculating the value of for a specific disease relies on secondary methods. There are many methods to calculate from mathematical models, few of which agree with each other and few of which produce the average number of secondary infections. Methods to calculate from theoretical models include the survival function, the next-generation method, the eigenvalues of the Jacobian matrix, the existence of the endemic equilibrium, and the constant term of the characteristic polynomial. can also be estimated from epidemiological data via the number of susceptibles at endemic equilibrium, the average age at infection, the final size equation and calculation from the intrinsic growth rate. For an overview, see Heffernan et al. [2].
Furthermore, there are many diseases that can persist with , while diseases with can die out, reducing the utility of the concept as a threshold. is also used as a measure of eradication for a disease that is endemic, but issues such as backward bifurcations, stochastic effects, and networks of spatial spread mean that an invasion threshold does not necessarily coincide with a persistence threshold. This results in a reduction of the usefulness of . For example, it is possible that a disease can persist in a population when already present but would not be strong enough to invade. Finally, the threshold value that is usually calculated is rarely the average number of secondary infections, diluting the usefulness of this concept even further.
In this paper, we outline the problems with and examine a number of alternatives that have been proposed. We include a worked example of malaria to demonstrate the many different results that the various methods give for the same model. Finally, we survey some of the recent uses of in the literature. The number of articles that use likely numbers in the tens of thousands, so an exhaustive review is not feasible. We have restricted ourselves to articles published since 2005 and which include interesting or novel explorations of .
2. Methods for Calculating
In this section, we identify some of the more popular methods (although by no means all) used to calculate . We also describe the limitations that each method presents and demonstrate one of the core problems with . Specifically, we address a key problem with : how do biologists make sense of it from mathematical models? (See, for example, the puzzled discussion in van den Bosch et al. [6].)
Although the ββ in is derived from βreproductiveβ, based on the original formulation of the concept as the average number of secondary infections, many thresholds have been denoted by ββ, even when they are not related to the average number of secondary infections. Thus, in keeping with the notation, we will use the notation to denote an -like surrogate associated with a particular method, symbolised by .
2.1. The Survival Function
The survival function is given by
The survival function has the advantage that it always produces the average number of secondary individuals infected by a single infected individual, in the same class. Thus, in Figure 1, where one human infects two mosquitoes, who each infect three humans, the survival function produces . This is the number of humans infected by a single infected human via mosquitoes; or, equivalently, the number of mosquitoes infected by a single mosquito via humans. The survival function is a generalised method of calculating the basic reproductive ratio that is not restricted to ODEs.
However, determining the individual probabilities can be cumbersome, especially if multiple states are involved. For a vector-borne infection such as malaria, with two infection states (human and mosquito), calculating the first probability involves determining the probability that a human infected at time 0 exists at time , the probability that a human infected for total time infects a mosquito and the probability that an infected mosquito lives to be age , where [2]. For diseases with more states, such as Guinea Worm disease, where there is a waterborne parasite, which can attach itself to copepods, which in turn are ingested by humans and which subsequently grow into an internal nematode, the calculations of the survival probabilities become unwieldy.
Thus, although this method always produces the correct , in practice, it is difficult to use. This is especially true for models with sufficient complexity, which are often those encountered most frequently.
2.2. The Jacobian
The Jacobian matrix is used to linearise a nonlinear system of differential equations. Around the disease-free equilibrium, the linear system will have the same stability properties as the nonlinear system if it is hyperbolic; that is, if no eigenvalues have zero real part. In particular, if all eigenvalues have negative real part, then the equilibrium is stable, whereas if there is an eigenvalue with positive real part, the equilibrium is unstable.
It follows that a threshold is , where is the largest eigenvalue of the Jacobian matrix (or the largest real part if the eigenvalues are complex). However, for a system of differential equations, this requires solving an th order polynomial, which may be impossible. Furthermore, rearranging the condition to produce a threshold is not a unique process and does not always produce the average number of secondary infections.
2.3. Constant Term of the Characteristic Polynomial
When , the constant term of the characteristic polynomial will be zero. However, the reverse is not true, as the polynomial could have both zero and positive roots. If the characteristic polynomial is then is a threshold if for all . More generally, the Routh-Hurwitz condition allows the coefficients to take on other signs, under certain restrictions, but the nonconstant coefficients all being positive is a sufficient condition. Another sufficient condition is that under the constraint (so that the largest eigenvalue at is 0).
This method is significantly easier to use than finding the largest eigenvalue, although verifying that necessarily corresponds to the largest eigenvalue can become complicated for some models. However, similar to the above, rearranging to produce is not a unique process and does not always produce the average number of secondary infections.
2.4. The Next-Generation Method
The next-generation method, developed by Diekmann et al. [8] and Diekmann and Heesterbeek [9], and popularised by van den Driessche and Watmough [10], is a generalisation of the Jacobian method. It is significantly easier to use than Jacobian-based methods, since it only requires the infection states (such as the exposed class, the infected class and the asymptomatically infected class) and ignores all other states (such as susceptible and recovered individuals). This keeps the size of the matrices relatively manageable.
However, the next-generation method suffers from a lack of uniqueness. In order to determine the matrices and (where accounts for the βnewβ infections and accounts for the transfer between infected compartments), biological insight must be used in order to decide which terms count as βnewβ infections and which terms are transfer terms. While this may seem intuitive for most models, it is easy to construct a counterexample Here, the term might represent, for example, new infections arising from vertical transmission, whereas the term might represent a disease-specific death rate. Although this construction is clearly arbitrary, it demonstrates that identifying βnewβ infections is not a unique process and relies on the modeller's judgement.
We would then have and thus This has the same threshold property as , but is clearly not the average number of secondary infections. In essence, the next-generation method is a mathematical generalisation of putting the negative values on one side and dividing (this is what is) so that the eigenvalue threshold at zero is transformed into an -like threshold at one. However, as we have seen, this does not produce a unique result.
van den Driessche and Watmough [10] note that other decompositions of and can be chosen, which lead to different values for . They claim that only one choice of is epidemiologically correct. However, this is not true, as the above example shows. Furthermore, it means that the definition of relies upon the judgement of the modeller as to what βepidemiologically correctβ means.
Furthermore, the next-generation method does not produce the number of humans infected by a single human if there is an intermediate host, but rather the geometric mean of the number of infections per generation.
For example, consider a mosquito-borne disease where humans infect two mosquitoes, while mosquitoes infect three humans, as shown in Figure 1. For convenience, label these and . Then the number of humans infected from a primary human (via mosquitoes) is . (This is also the value calculated by the survival function.)
However, the next-generation method would calculate , which is a weighted average () of the number of infectives each individual produces in the next infection event. While mathematically sound, it is questionable whether this is biologically meaningful.
This example could be extended. Consider a three-stage disease, such as tularemia, where ticks may transmit between humans and livestock, but humans may also be infected by eating livestock. Suppose that a single human directly infects two ticks (so ), each tick infects four animals (so ) and each animal infects three humans (so ). Then, a single human has resulted in infected humans. However, , since there are three infection stages. As the number of infection stages increase, becomes a progressively higher surd.
The next-generation method is likely the most frequently used method to calculate . It has been used extensively to calculate -like values from host-vector models (see Wonham et al. [11], Gaff et al. [12] or Gubbins et al. [13]). However, it does not produce the number of newly infected individuals in the same infection class and does not always produce the average number of secondary infections.
2.5. The Graph-Theoretic Method
In de Camino-Beck et al. [14], a graph-theoretic method for calculating is given. Starting from the definition of , they derived a series of rules for reducing the digraph associated with to a digraph with zero weight, from which . The rules are as follows.
Rule. To reduce the loop to β1 at node , every arc entering has weight divided by .
Rule. For a trivial node on a path , the two arcs are replaced by with weight equal to the product of weights on arc and . Weights on multiple arcs are added.
The graph-theoretic method has the advantage that it avoids calculation of , which may be complicated for large systems. However, it always produces the same threshold value as the next-generation method.
They considered the simple vector-host model where , , , and are proportions of infective hosts, infective vectors, susceptible hosts, and susceptible vectors, respectively, is the host birth and death rate, is the host recovery rate, is the vector birth and death rate, and and are disease transmission coefficients. The disease-free equilibrium is . They found matrices From this, they produced a digraph of and concluded that
However, this method still contains the same issues as the next-generation method. The value calculated is not the number of humans infected by a single human, but rather the (less biologically meaningful) geometric mean of the number of humans infected by the vector and the number of vectors infected by a human. Furthermore, the creation of and is not a unique process, as we have seen above, and does not always produce the average number of secondary infections.
For example, consider the mathematically equivalent model with matrices
The graph reduction then proceeds as in Figure 2. It follows (either by the graph-reduction method or using the conventional next-generation method) that Thus, the graph-theoretic method inherits the existing problems from the next-generation method.
2.6. Existence of the Endemic Equilibrium
can also be calculated from the endemic equilibrium, in the case where there is a bifurcation at such that the endemic equilibrium does not exist for . The existence of the endemic equilibrium is, thus, a threshold that can be rearranged to produce an -like surrogate.
However, for many models, calculating the endemic equilibrium can be quite difficult. Furthermore, in the case of a backward bifurcation, the endemic equilibrium still exists for (in fact, two endemic equilibria exist in this case, one of which is stable and the other unstable). It follows that the endemic equilibrium is not a useful general method for calculating and does not always produce the average number of secondary infections.
van den Bosch et al. [6] use the endemic equilibrium to determine a threshold given by so that if (where represents infected plants at the endemic equilibrium, is the transmissibility, is the death rate, is the fraction of infectious seeds, and is the total plant population density). They show that two different rearrangements of this inequality give and note that either would suffice as an value, but were unable to resolve the question of which was appropriate.
2.7. Summary of Methods
In summary, there are many methods available for calculating , but few of them agree with each other and almost none reliably calculate the average number of secondary infections in a wholly susceptible population. The only method which does is the survival function, but this method is too cumbersome to be used except for the simplest of models.
3. A Worked Example
In this section, we take a sample model and apply the various methods for calculating to it. We show that each method can produce a different result, all of which have the property that they have an outbreak threshold at but otherwise bear little relation to one another.
Consider the following model for malaria: Humans may be susceptible or infected ( and , resp.), while mosquitoes may be susceptible or infected ( and , resp.). The birth rates are for humans and for mosquitoes. The background death rates are and for humans and mosquitoes, respectively, while the disease-specific death rate for humans is . The transmission rate from mosquitoes to humans is , while the transmission rate from humans to mosquitoes is . See Figure 3.
The Jacobian matrix at the disease-free equilibrium is where and are the equilibrium values of susceptible humans and mosquitoes, respectively. Thus, the nontrivial part of the characteristic polynomial satisfies
Using the Jacobian method, the largest eigenvalue is Rearranging , we have
Conversely, since the nonconstant coefficients of the characteristic polynomial are positive, all eigenvalues will be negative if the constant term of the characteristic polynomial is positive, whereas there will be an eigenvalue with positive real part if the constant term is negative. Rearranging, we have
This is not the only rearrangement possible, so, like the nonuniqueness of the next-generation method, it is possible to rearrange by adding and subtracting arbitrary constants. Thus, for example, is also a measure of disease spread with the property that the disease persists if and is eradicated if .
Other rearrangements are possible, so long as the threshold at 0 is transformed into a threshold at 1. Thus, has the same threshold property. A generalised version of this formulation was used in Smith? et al. [15] called .
The endemic equilibrium satisfies It follows that there will be a biologically meaningful endemic equilibrium if , or since and . As above, this is not the only rearrangement possible.
The next-generation matrices are since the next-generation method only considered the two infected classes. Then, using the next-generation method,
Thus, for the same model, , , , , and are all distinct. Although in this case, this does not hold in general. Note that and produce the number of infected humans per infected human (or, equivalently, the number of infected mosquitoes per infected mosquito), but this is also not necessarily true in general. Figure 4 illustrates how these -like expressions vary with , when all other parameters are held constant.
Anecdotally, we have found that most biologists believe that there is one and only one value for a given disease and, in order to verify an is correct, all that is required is to check that it has a threshold at 1. We have demonstrated here that is not a unique concept. Indeed, it is straightforward to construct simple variations that satisfy the threshold criterion at 1. Furthermore, since multiple methods calculate multiple thresholds, it follows that the vast majority of them cannot be the average number of secondary infections. However, the nonuniqueness is not the only problem associated with the concept, as we will demonstrate.
4. Backward Bifurcations
Backward bifurcations occur when multiple stable equilibria coexist for . This presents a serious complication when a disease is already endemic, since lowering the basic reproduction number below 1 may no longer be a viable control measure; hence, different prevention and control measures may have to be considered. In particular, a backward bifurcation makes the system more complicated, since the behaviour now depends on the initial conditions.
A backward bifurcation at may result in persistence of the disease when . In this case, the disease will always persist for . However, there is a point such that the endemic equilibrium exists for and a third, unstable, equilibrium also exists between the two stable equilibria. Hence, an endemic disease is only eradicated if . For , the outcome depends on initial conditions. If the disease is still in its early stagesβthat is, if the initial conditions are sufficiently smallβthen the system will approach the disease-free equilibrium and the disease will be eradicated. However, if the initial conditions are large, then the system will approach the endemic equilibrium and the disease will persist. Thus, a backward bifurcation prevents the system from switching to the disease-free equilibrium as soon as is reached. See Figure 5(a).
(a)
(b)
Several mechanisms have been shown to lead to backward bifurcation in epidemic models, but backward bifurcations in compartmental models have only recently attracted serious research attention.
GΓ³mez-Acevedo and Li [16] investigated a mathematical model for human T-cell lymphotropic virus type I (HTLV-I) infection of CD4+ T cells that incorporates both horizontal transmission through cell-to-cell contact and vertical transmission through mitotic division of infected T cells. They assumed that a fraction of the infected cells survive the immune system attack after the error-prone viral replication. Under the biologically sound assumptions that the fraction should be very low and the rate of the mitotic division should be high, their model has a bifurcation that predicts persistent infection for an extended range of the basic reproduction , where . This model undergoes a backward bifurcation as increases: multiple stable equilibria exist for an open set of parameter values, when the basic reproduction number is below one.
Safan et al. [17] studied an epidemiological model under the assumption that the susceptibility after a primary infection is times the susceptibility before a primary infection. They present a method for determining the control effort required to eliminate an infection from a host population when subcritical persistence may occur. This effort can be interpreted as a reproduction number but is not necessarily the basic reproduction number.
For , this model exhibits backward bifurcations, where is the death rate and is the recovery rate. For such models, the authors presented a method for determining the minimum effort required to eradicate the infection from the endemic steady state if one concentrates on control measures affecting the transmission rate constant.
Garba et al. [18] presented a deterministic model for the transmission dynamics of a single strain of dengue by realistically adopting a standard incidence formulation and allowing dengue transmission by exposed humans and vectors. The model was extended to include an imperfect vaccine for dengue. A backward bifurcation was observed in both models. This makes no longer sufficient for effectively controlling dengue in a community. However, this phenomenon can be removed by replacing the standard incidence function in the model with a mass-action formulation.
Reluga et al. [19] proposed a series of epidemic models for waning immunity that can be applied in many different settings. With biologically realistic hypotheses, they found that immunity alone never creates a backward bifurcation. However, this does not rule out the possibility of multiple stable equilibria, which can be shown by a counterexample of the forward bifurcation at . See Figure 5(b).
5. in Spatial Contexts
When is considered in a spatial context, many of its properties fail to hold. In particular, diseases with can fail to persist, depending on the nature of the spatial transmission. Since many diseases are spatially dependent, this further limits the utility of . As increases beyond 1, the probability of disease invading the initially infected host group increases, but additional criteria are important to determining the probability of the spreading of the disease to other groups.
5.1. Networks
Green et al. [20] related deterministic mean-field models to network models, taking into account the manner in which the contact rate and infectiousness change over time. For a node with exactly connections, the expected number of secondary infections at infection age is given by where is the transmissibility and is the mean number of connections per node. This model applies when the rate of infectious contact is independent of .
If there is a constant rate of generation of new cases, then the expected number of secondary infections in an infectious period of length is given implicitly by where denotes the contact per susceptible neighbour and denotes the infectiousness at time after infection.
Kao [21] showed that novel pathogens may evolve towards a lower , even if this results in pathogen extinction. This is because the presence of exploitable heterogeneities, such as high variance in the number of potentially infectious contacts, increases ; thus, pathogens that can exploit heterogeneities in the contact structure have an advantage over those that do not. The exploitation of heterogeneities results in a more rapid depletion of the potentially susceptible neighbourhood for an infected host. While the low strategy is never evolutionarily stable, invading strains with higher will also converge to the low strategy if not sufficiently different from the resident strain. This is in contrast to the conventional belief that the emergence of novel pathogens is driven by maximisation of .
In a randomly mixed epidemiological network, can be approximated by where and are, respectively, the number of inward and outward βtruly infectiousβ links per node; the angled brackets represent the expected value of the relevant quantity [22].
Meyers [23] showed that in a contact network framework, where is the mean probability of transmission between individuals and and are the mean degree and mean square degree of the network. Here, depends explicitly on the structure of the network (i.e., on and ). A single pathogen may, therefore, have very different transmission dynamics depending on the population through which it spreads. If two networks have the same mean degree, , then the one with the larger variance in degree, , will be more vulnerable to the spread of disease.
In compartment models, infected hosts are assumed to have potentially disease-causing contacts with random individuals from the population according to a Poisson process that yields an average contact rate of per unit time. The mass-action assumption of compartmental models is tantamount to assuming that the underlying contact patterns form a random graph with a Poisson degree distribution. Estimates of that assume a mass-action model may, therefore, be invalid for populations with non-Poisson contact patterns and, in particular, will underestimate the actual growth rate of the disease in highly heterogeneous networks.
5.2. Individual-Level Models
Rahmandad and Sterman [24] noted that if in an agent-based model then, due to the stochastic nature of interactions, it is possible that no epidemic occurs or that it ends early if, by chance, the few initially contagious individuals recover before generating new cases.
Schimit and Monteiro [25] showed that in an individual-level model, cannot be uniquely determined from some features of transient behaviour of the infective group. The value of can be unambiguously determined from the asymptotical stable stationary concentrations, but this relies on waiting for the system to reach its permanent regime, which is not feasible in practice. The same value of can be associated to networks with distinct values of clustering coefficients and average shortest path length. This result can affect the evaluation of the effectiveness concerning different strategies employed for controlling a disease. Because distinct values of topological properties can produce the same value of in a model considering the spatial structure of the contact network, it is difficult to evaluate the effective contribution of each control measure. This is because the correspondences among and the topological properties of the contact network are not one-to-one.
5.3. Metapopulation Models
Cross et al. [26] showed that when is based on data collected from simulated epidemics mimicking epidemiological contact-tracing data, can be substantially greater than one and yet not cause a pandemic. In populations with social or spatial structure, a chronic disease is more likely to invade than an acute disease with the same , because it persists longer within each group and allows for more host movement between groups.
Under the settings where the rate of host population turnover was negligible relative to the rate of disease processes of infection and recovery, they showed that was insufficient for disease invasion when the product of the average group size and the expected number of between-group movements made by each individual while infectious was less than 1.
Smith? et al. [15] examined a metapopulation model with travel between two regions, with reproductive ratios and for each region in the absence of travel, and and when only susceptibles travel, but infectives do not. They showed that if both and , then there are conditions on the travel of susceptible such that and . Thus, a disease which would otherwise be eradicated in both regions could be sustained in one of the regions if there were sufficient travel of the susceptibles (not infectives). Furthermore, if and , then there are conditions on the travel of susceptibles such that and . Thus, if one region sustains the disease on its own, while the other does not, then sufficient travel of susceptibles (not infectives) could sustain the disease in both regions.
5.4. Partial Differential Equation Models
Althaus et al. [27] examined an age-dependent partial differential equation model of in-host HIV infection. They showed that where the denominator is the Laplace transform of the generation time distribution and is the growth rate. They found that estimates for were generally smaller than those derived from the standard model when the generation time was taken into account.
6. Stochastic Effects
When stochastic effects, such as those inevitably found in nature, are included, the threshold at may be disturbed. This includes assumptions about the distribution of transition times (assumed to be exponential in most models) as well as variations in individual parameters.
Heffernan and Wahl [28] derived improved estimates of for situations in which information about the dispersal of transition times is available to the clinical or epidemiological practitioner. Rather than rederiving for a number of models (SIR, SEIR, etc.), they introduce a βcorrection factorβ, : the ratio of when the lifetimes are nonexponentially distributed to the value of that would be calculated assuming exponential lifetimes. They were able to derive limiting values of and used this to gauge the sensitivity of to dispersion in the underlying distributions.
By combining the movement of hosts, transmission with groups, recovery from infection and the recruitment of new susceptibles, Cross et al. [29] expanded the earlier analysis of Cross et al. [26] to a much more broader set of disease-host relationships, exploring settings where the duration of immunity ranges from transient to lifelong or where the demographic processes occur on comparable (or faster) timescales to disease processes. The focus of this study was to investigate how recruitment of susceptibles affects disease invasion and how population structure can affect the frequency of superspreading events (SSEs). They found that the frequency of SSEs may decrease with the reduced movement and the group sizes due to the limited number of susceptibles available.
The hierarchical nature of disease invasion in host metapopulations is illustrated by the classification tree analysis of the model results. Firstly, the pathogen must effectively transit within a group (), and then, the pathogen must persist within a group long enough to allow for movement between different groups. Hence, the infectious period, group size and recruitment of new susceptibles are as important as the local transmission rates in predicting the spread of pathogens across a metapopulation. It should be noted that in 35% of simulations when was greater than one, the disease failed to invade.
Tildesley and Keeling [30] examined whether was a good predictor of the 2001βUK Foot and Mouth disease. They concluded that explained just 29.3% of the standard deviation of the epidemic impact. They also noted that did not act as a threshold: at the value of , only 20% of initial seedings generated epidemics; this probability increased to around 50% for the largest reasonable values. When heterogeneities exist in the population, infection is most likely to become focused within the high-risk individuals who are both more susceptible and more infectious. This highlights the stochastic nature of the disease in its early stages and the dependence of the ensuing epidemic on favourable local conditions in the neighbourhood of the initial infection.
The probability of extinction, assuming exponentially distributed infectious periods, was when an epidemic began with a single infected individual. Thus, if , then extinction occurs with probability 1, but if then extinction occurs with some probability. See Chiang [31, Chapter 4] for more discussion.
7. Failures
In this section, we note a variety of problems with that are not covered in the previous sections. These include problems with the underlying structure of compartment models, the mismatch between an individual-based parameter and a population-level compartment model, and the failure of to accurately measure an outbreak of a new disease.
Breban et al. [32] argued that in order to associate an to a model of ODEs, an individual-level model (ILM) which is compatible to the ODE model must be developed; only then can the of the ILM be unambiguously calculated. These ILMs are growing (not static) network models, with individuals added to a network of who infected whom based on global or local network rules. Then, is computed as the limit of the average number of outgoing links of individuals in a node that no longer accepts new links, as time goes to infinity. They showed that a broad range of values were compatible with a given ODE model.
For example, consider the basic model where is the transmissibility and is the disease death rate; note that the transmission term is equivalent to the standard term when the depletion of susceptibles is negligible so that .
The expected value from this model using other methods (such as the next-generation method) is The corresponding ILM consists of an infection rule, where an individual joining the infectious pool is infected by an infectious individual who is uniformly randomly selected, and a removal rule, where a uniformly randomly selected individual leaves the infectious pool. The flow of newly infected individuals is . Thus, the flow per already-infected individual is . Since the removed individuals are randomly sampled from the infectious individuals, the average length of the infectious period equals the time expectation of the infectious period, which is (rather than , as calculated from the next-generation method). It follows that , which is independent of the transmission and death rates from the corresponding ODE. Thus, in this case, does not signal epidemic growth as anticipated from other methods.
Roberts [3] noted three fundamental properties commonly attributed to : (i) that an endemic infection can persist only if , (ii) provides a direct measure of the control effort required to eliminate the infection, and (iii) pathogens evolve to maximise their value. He demonstrated that all three statements can be false. The first, as we have noted, can fail due to the presence of backward bifurcations. The second can fail when control efforts are applied unevenly across different host types (such as a high-risk and a low-risk group), since is determined by averaging over all host types and does not directly determine the control effort required to eliminate infection.
The third can fail when two pathogens coexist at a steady state that exists and is stable whenever both single-pathogen steady states exist but are unstable. In this case, the order in which the pathogens are established in the host population matters. The established parasite has a role in determining a modified carrying capacity and the pathogen with the largest basic reproductive ratio does not necessarily exclude the other.
Breban et al. [33] showed that two individual-level models having exactly the same expectations of the corresponding population-level variables may yield different values. They showed that obtaining from empirical contact-tracing data collected by epidemiologists and using this as a threshold parameter for a population-level model could produce misleading estimates of the infectiousness of the pathogen, the severity of an outbreak, and the strength of the medical and/or behavioural interventions necessary for control. Thus, measuring through contact tracing (as generally occurs during an outbreak investigation) may not be a useful measure for determining the strength of the necessary control interventions.
Many different individual-level processes can generate the same incidence and prevalence patterns. Thus, assigning a meaningful value to an ODE model without knowledge of the underlying disease transmission network may be impossible. Only an epidemic threshold parameter can be used to design control strategies. Since fails to possess this threshold quality, its usefulness may be vastly overstated.
Meyers [23] compared the theoretical calculation of with observed SARS data from China and showed that estimates of seemed incompatible. The basic reproductive rate has two critical inputs: (1)intrinsic properties of the pathogen that determine the transmission efficiency per contact and the duration of the infectious period,(2)the patterns of contacts between infected and susceptible hosts in the population.
While the first factor may be fairly uniform across outbreaks, the second may depend significantly on context, varying both within and among populations. The problem with the SARS estimates stems from the mass-action assumption of compartmental models; that is, that all susceptible individuals are equally likely to become infected. When this assumption does not hold, the models may yield inaccurate estimates or estimates that do not apply to all populations. estimates for SARS in the field were based largely on outbreak data from a hospital and a crowded apartment building, with anomalously high rates of close contacts among individuals.
The author suggested that it might be inappropriate to extrapolate estimates for from specific settings such as these to the population at large. Conversely, since the population at large is also unlikely to satisfy mass-action requirements, it may also be concluded that is not a meaningful estimate of disease spread.
8. Alternatives to
A number of alternatives have been offered over the years, due to recognised problems with . Older examples include the critical community size [5], the group-level reproductive number, [34], the type reproduction number [35] and the basic depression ratio [36]. See Heffernan et al. [2] for a summary of these methods. In this section, we review some of the more recent alternatives to .
The actual reproduction number, , is defined as a product of the mean duration of infectiousness and the ratio of incidence to prevalence [37β39]. coincides with when the transmission probability is constant but accounts for the more general situation when the transmission probability varies as a function of infection age (as happens in diseases such as HIV/AIDS).
The effective reproduction number measures the number of secondary cases generated by an infectious case once an epidemic is underway. In the absence of control measures, , where is the proportion of the population susceptible [40, 41]. The estimation of reproductive numbers is typically an indirect process because some of the parameters on which these numbers depend are difficult or impossible to quantify directly. The effective reproductive number satisfies , with equality only when the entire population is susceptible [40].
The effective reproduction number is of practical interest, since it is time dependent and can account for the degree of cross-immunity from earlier outbreaks. However, since it is based on the basic reproduction number, the effective reproduction number inherits many of the issues from .
Breban et al. [32] proposed , the average number of secondary infections over the infectious population. The average number of secondary infections of actively infectious individuals, , is computed as the average number of outgoing links of a node in the infected compartment at time . Then, when the limits exist. Unfortunately, is never defined in the paper, limiting the usefulness of this formulation of . However, by analogy with , is the average number of outgoing links of a removed compartment at time (Romulus Breban, personal communication).
For the SI model (32), under the assumption that every infection is uniquely assigned as a secondary infection for either a removed or an infected individual, where is the cumulative number of infected individuals that occur in the time interval and is the cumulative number of removed individuals.
Their definition of evaluates the average number of secondary cases over removed individuals as the distribution of secondary cases becomes stationary. This definition does not imply a particular individual-level model; it depends exclusively on the structure of the disease-transmission network [42]. However, is not measured at the start of an epidemic, which may limit its usefulness during an initial outbreak.
Grassly and Fraser [43] demonstrated that standard epidemiological theory and concepts such as do not apply when infectious diseases are affected by seasonal changes. They instead define where is the transmission parameter at time and is the average duration of infection. Thus, can be interpreted as the average number of secondary cases arising from the introduction of a single infective into a wholly susceptible population at a random time of the year.
They noted that the condition is not sufficient to prevent an outbreak, since chains of transmission can be established during high-infectious seasons if . However, is both necessary and sufficient for long-term disease extinction.
Meyers [23] noted that in the case of networks, estimating the average transmissibility may be more valuable than . This means reporting not only the number of new infections per case, but also the total estimated number of contacts during the infectious period of that case. Given the primary role of contact tracing in infectious disease control, these data are often collected. Unlike , can be justifiably extrapolated from one location to another even if the contact patterns are quite disparate.
Instead of , the author offered a number of alternatives to determine whether an outbreak will occur, based on network modelling. These include the probability that an individual will spark an epidemic and the probability that a disease cluster will spark an epidemic.
The probability that an individual with degree will spark an epidemic is equal to the probability that transmission along at least one of the edges emanating from that vertex will lead to an epidemic. For any of its edges, the probability that the disease does not get transmitted along the edge is . The probability that disease is transmitted to the attached vertex but does not proceed into a full-blown epidemic is , where is the probability that a secondary infection does not spark an epidemic. Thus, the probability that an individual will spark an epidemic is .
The probability that an outbreak of size sparks an epidemic is where is the relative frequency of vertices of degree in the network.
Kao et al. [22] defined an epidemiological network contact matrix whose elements are either 1 or 0, depending on whether an infectious contact between nodes and is possible. The spectral radius of is an alternative approximation for , which can be calculated via a weighted version of (28).
This explicitly accounts for the full contact structure of the network, but the evaluation of extremely large, reasonably dense matrices (some highly active nodes may have hundreds of potentially infectious links) is difficult and time consuming, particularly when this evaluation process must be repeated multiple times. However, comparisons between the two approximations for subsets of a sheep network with several thousand nodes show little difference between and the spectral radius of (typically less than 5%).
Kamgang and Sallet [44] used the special structure of Metzler matrices (real, square matrices with nonnegative off-diagonal entries) to define , an analytical threshold condition. is a function of the parameters of the system such that if , the disease-free equilibrium is locally asymptotically stable, and if , the disease-free equilibrium is unstable. has an association with , although it is a stronger condition; however, it has no direct biological interpretation. The algorithm for deriving , although highly mathematical in nature, allows computation of a threshold for high-dimensional epidemic models.
Huang [45] defined four reproductive numbers associated with four types of transmission patterns, each depending on , the ratio of the mean infectious period to the mean latent period. These four reproduction numbers are the following: (1), the minimal reproductive number associated with the slowest latency process and the fastest recovery process,(2), the middle reproductive number associated with mean latency and recovery processes,(3), the maximal reproductive number associated with the fastest latency process and the slow recovery process,(4), the largest reproductive number associated with the fastest latency process and the extremely slow recovery process.
All four reproduction numbers are strictly increasing functions of and satisfy These numbers allow a disease to be classified as mild () or severe ().
Hosack et al. [46] noted that does not necessarily address the dynamics of epidemics in a model that has an endemic equilibrium. They used the concept of reactivity to derive a threshold index for epidemicity, , which gives the maximum number of new infections produced by an infective individual at a disease-free equilibrium. They also showed that the relative influence of parameters on and may differ and lead to different strategies for control.
If is derived from the next-generation operator so that then the threshold for endemicity is such that the system is reactive when and nonreactive when . describes the transitory behaviour of disease following a temporary perturbation in prevalence. If the threshold for epidemicity is surpassed, then disease prevalence can increase further even when the disease is not endemic. This suggests that epidemics can occur even in areas where long-term transmission cannot be maintained.
Reluga et al. [47] defined 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. It is calculated as where is the discount rate, is the identity matrix, and and are from the next-generation matrix decomposition. combines properties of both the basic reproductive number and the ultimate proliferation rate although it also inherits the nonuniqueness problems from the next-generation method.
Nishiura [4] developed a likelihood-based method for estimating without assuming exponential growth of cases and offers a corrected value of the actual reproduction number. The author noted that is extremely sensitive to dispersal of the progression of a disease or variations in the underlying epidemiological assumptions.
9. Discussion
Despite being a crucial concept in disease modelling, with a long history and frequent application, is deeply flawed. It is not a measure of the number of secondary infections, it is never calculated consistently, and it fails to satisfy the threshold property. Rarely has an idea so erroneous enjoyed such popular appeal. Why, then, are we so attached to the concept?
The answer is that , despite all its flaws, is all that we have. No other concept has so effectively transcended mathematics, biology, epidemiology, and immunology. No other concept is so general that it can be understood in terms of compartment models, network models, partial differential equations, and metapopulation models. βThe number of secondary infectionsβ has an intuitive appeal that outlasts even the inaccuracy of that statement when applied to the concept.
The threshold nature of is used to monitor and control severe real-time epidemics; control measures are often concluded if [48], making the problems with more than just theoretical. Due to the inconsistencies in calculation, different diseases cannot be compared unless the same method was used to calculate ; if HIV has an of 3 and swine flu has an of 4, we cannot conclude that swine flu is worse than HIV if different methods were used to determine these values. All we can conclude is that both diseases have ; however, as we have demonstrated, this does not necessarily guarantee disease persistence.
Of the many different methods used to calculate , only the survival function reliably calculates the average number of secondary infections; however, this method is too cumbersome to use in most practical settings. The next-generation method is probably the most popular, yet it suffers from uniqueness problems and does not cope well with more than one disease state. Since is rarely measured in the field, it instead relies upon after-the-fact calculations to determine the strength of disease spread [2]. This limits its usefulness even further.
Policy decisions are being based upon the concept, with limited understanding of the complexity and errors that exist in the very structure of the concept. Funding decisions about where money should be best spent are based on estimates of , resources are directed towards one disease over another, and monitoring programmes are abandoned, their objectives only half-realised, because of . Lives may be saved or lost, based on this imperfect and inconsistent measure.
is a quantity that relates to the initial phase of an epidemic. This makes practical sense in terms of disease prevention. However, it is also used to guide eradication efforts when a disease is endemic. Some methods derive an eradication threshold from an equilibrium value that may not be attained for a very long time. This suggests that different measures are needed in different stages of an epidemic in order to characterise transmissibility and to guide intervention strategies. When a new pathogen emerges, a quantity describing the initial spread is useful. When a disease is endemic, a quantity that applies to the long-term equilibrium is more appropriate.
We note that although the definition of is broad, it is not a universal quantity that applies to all settings. In different settings, one should use a quantity that satisfies the following properties: that an endemic equilibrium only persist if and that provides a direct measure of the control effort required for eradication. The contact structure of the population, the variation of risk factors and the order of establishing parasites (if applicable) should accompany the identification of a meaningful .
What is urgently needed is a simple, but accurate, measure of disease spread that has a consistent threshold property and which can be understood by nonmathematicians. If is to be used, it must be accompanied by a declaration of which method was used, which assumptions are underlying the model (e.g. mass-action transmission) and evidence that it is actually a threshold, with no backward bifurcation. Without such caveats, the concept of will continue to fail.
Appendix
A. for Specific Diseases
In this appendix, we illustrate some recent calculations of for various diseases from the past few years, noting in particular specific values for that have been given. This illustrates the wide variety of values that are presented as being βtheβ value for a specific disease.
A.1. HIV
BacaΓ«r et al. [49] developed a mathematical model of the HIV/AIDS epidemic in Kunming, China. They considered needle sharing between injecting drug users and commercial sex between female sex workers and clients. The basic reproductive ratio is expressed as where the superscript βIDUβ represents injecting drug users and βsexβ represents sex workers. They estimate and .
Abu-Raddad et al. [50] examined the population-level implications of imperfect HIV vaccines. They divided the population into two groups: unvaccinated and vaccinated, both of whom could become infected (albeit at different rates, depending on the efficacy of the vaccine). By forming two independent next-generation matrices, they derived and , the basic reproductive ratios for the unvaccinated and vaccinated populations, respectively. They note that the condition may not be sufficient for disease eradication, due to the possibility of a backward bifurcation. The two formulas are equivalent in the absence of the vaccine, assuming that risky behaviour does not change. They estimate in the absence of vaccine protection.
Gran et al. [39] argued that was inappropriate for long-lasting infections and instead defined the actual reproduction number . This is a time-dependent measure defined as the average number of secondary cases per case to which the infection was actually transmitted during the infectious period in a population. Thus, where is the incidence rate at time , is the prevalence at time , and is the average length of the infectious period. They illustrated their results for a staged HIV/AIDS progression model for homosexual men in the UK and showed that was 13.81 and 17.01 in 1982 and 1983, respectively, but was much closer to 1 by the mid 1990s, due to decreases in risk behaviour.
Smith? et al. [15] used the nonuniqueness of -like thresholds to evaluate the feasibility of eradicating AIDS using available resources. They developed a linear metapopulation model that has the same eradication threshold as more realistic models. They defined the threshold as where is the maximal real part of all eigenvalues of , the Jacobian matrix evaluated at the disease-free equilibrium; represents the change of status of infected individuals in different regions. The model has the property that the disease will be eradicated if but will persist if . The paper stresses that this quantity should be understood only as a threshold condition for eradication; the model from which it derives does not quantify the transient dynamics, the timecourse of the infection or the prevalence of the disease. However, they showed that the simpler version of the model has the same eradication threshold () as more realistic, nonlinear models. The mathematical analysis results, together with a simple cost analysis, used the threshold nature of to show that eradication of AIDS is feasible, using the tools that we have currently to hand, but action needs to occur immediately and globally.
A.2. Influenza
Chowell et al. [51] used four methods to calculate using data from the 1918 influenza pandemic in California. The four methods involved (1) an early exponential growth rate, (2) an SEIR model, (3) a more complicated SEIR model with asymptomatic and hospitalised cases, and (4) a stochastic SIR model with Bayesian estimation. These yielded average values of (1) 2.98, (2) 2.38, (3) 2.2, and (4) 2.1 during the first 17 days of the epidemic and then 2.36 for the entire autumn wave. In particular, uncertainty during the early part of the epidemic leads to a wide range of estimates for (0.5β3.5, 95% confidence interval), but uncertainty decreases with more observations.
Wallinga and Lipsitch [52] investigated how generation intervals shape the relationship between growth rates and reproductive numbers. For new emerging infectious diseases, the observed exponential epidemic growth rate can indirectly infer the value of the reproductive number. The existence of a few different equations that relate the reproductive number to the growth rate makes such inference ambiguous. It is unclear which of these equations might apply to new infection. The authors showed that these different equations differ only with respect to their assumed shape of the generation-interval distribution. They found that the shape of the generation-interval distribution determines which equation is appropriate for inferring the reproductive number from the observed growth rate. By assuming all generation intervals to be equal to the mean, they obtained an upper bound on the range of possible values that the reproductive number may attain for a given growth rate. They also showed that it is possible to obtain an empirical estimate of the reproductive number by taking the generation-interval distribution equal to the observed distribution.
The authors illustrated the impact that various assumptions about the shape of the generation interval distribution may have on the estimated value of the reproduction numbers for a given growth rate using human infections with influenza A virus. Observed generation intervals for influenza A in a Japanese household study excluding possible coprimary and tertiary cases [53] estimated a mean generation interval of 2.85 days and a standard deviation of 0.93 days. Without specific assumptions about the shape of the generation-interval distribution, they found that the reproductive number of influenza A is larger than , but smaller than .
Chowell and Nishiura [54] reviewed quantitative methods to estimate the basic reproductive ratio of pandemic influenza. They use the intrinsic growth rate to estimate in the range from 1.36 to 2.07. They defined the effective reproduction rate where is the reproductive power at time and infection age at which an infected individual generates secondary cases. If contact and recovery rates do not vary with time, then this simplifies to where is the population of susceptibles at time .
They also outlined statistical methods of estimation and showed that the expected number of secondary transmissions is given by where the pattern of secondary transmission follows a discrete probability distribution with secondary transmissions. This approach accounts for demographic stochasticity, which can be crucial during the initial phase of an epidemic when the number of infected individuals is small. This has been applied to observed outbreak data for avian influenza, where extinction was observed before growing to a major epidemic [55].
Tiensin et al. [56] used flock-level mortality data and statistical back-calculation methods to estimate from an SIR model for the 2004 avian influenza epidemic in Thailand. They calculated to be between 2.26 and 2.64, depending on the length of the infectious period.
A.3. Malaria
Smith et al. [57] revisited the basic reproduction number of malaria and its implications for malaria control. Their estimates of are based on two more commonly measured indices called the entomological inoculation rate, which is the average number of infectious bites received by a person in a year, and the parasite ratio, which is the prevalence of malaria infection in humans. They made 121 estimates of for Plasmodium falciparum malaria in African populations. The estimates (which do not come from a mathematical model) range from around one to more than 3000; the median is 115 and the interquartile range is 30β815.
They then revised the formula of to adopt potential sources of bias when biting is heterogeneous and when there is some transmission-blocking immunity. Under the assumption that the transmission-blocking immunity develops, the estimates of ranged from below one to nearly 11,000, with a median of 86 and an interquartile range of 15β1,000. This is because transmission-blocking immunity and heterogeneous biting skew the probability that a mosquito becomes infected, per bite. They also showed that in small human populations, the classical formulas of approximate transmission when counting infections from mosquito to mosquito after one parasite generation, but they overestimate it from human to human.
A.4. West Nile Virus
Wonham et al. [11] examined the conflicting outbreak predictions for West Nile virus, showing how the choice of transmission term affects . They concluded that some transmission terms apply biologically only at certain population densities and showed that six common North American bird species would be effective outbreak hosts. All seven models considered shared a similar compartment structure, but with varying factors, such as incubation periods, loss of immunity, age structure, and so forth. They created a core model synthesised from all seven models and used the next-generation method to calculate where is the transmissibility of the reservoir, is the death rate of the reservoir, is the death rate of the vector, is the total vector population at the disease-free equilibrium and is the total reservoir population at the disease-free equilibrium.
The first term under the square root sign represents the number of secondary reservoir infections caused by one infected vector. The second term under the square root sign represents the number of secondary vector infections caused by one infected reservoir. The square root gives the geometric mean of the two terms (see Section 2). Using reported minimum, mean and maximum estimates, 10,000 Monte Carlo estimates were run to generate an estimated distribution of . Mean numerical values for for house sparrows ranged from 25β1200.
However, these values are enormous, implying that a single sparrow would infect up to 1.44 million sparrows (i.e., ) during its infectious lifetime. Cruz-Pacheco et al. [58], for example, estimated for the house sparrow, using data from Komar et al. [59].
A.5. Cholera
Hartley et al. [60] developed a compartment model of cholera incorporating both a hyperinfective state and a lower infective state. They used the average age of infection to estimate from data for four cholera outbreaks in developing countries during the nineties. values ranged from 3.1 in Indonesia (1993β1999) to 15.3 in Pakistan (1990β1995), with an average of 8.7. They then used the next-generation method to develop a theoretical estimate for that incorporates hyperinfectivity. They show that this raises from approximately 3.2 when there is no hyperinfectivity to 18.2 when the contact with hyperinfected water is similar to contact with nonhyperinfected water, but that could be significantly larger with more frequent contact with hyperinfected water.
A.6. Dengue
Nishiura [61] clarified the contributions of mathematical and statistical approaches to dengue epidemiology. This highlighted the practical importance of the basic reproduction number, , in relation to the critical proportion of vaccination required to eradicate the disease in the future. The author illustrated three different methods to estimate : (i) the final size equation, (ii) the intrinsic growth rate, and (iii) age distribution, with published estimates and examples. The author pointed out that, although the estimates of most likely depend on the ecological characteristics of the vector population, it would be appropriate to assume a serotype-nonspecific estimate for is approximately 10, at least when planning vaccination strategies in endemic areas.
A.7. Rift Valley Fever
Gaff et al. [12] modelled Rift Valley fever, a mosquito-borne disease affecting both humans and livestock that currently exists in developing countries but has the potential to be introduced to the western world. The disease can be transmitted both horizontally (from host species to host species via mosquitoes) or vertically (via parent-to-offspring transmission in mosquitoes). They expressed the basic reproductive ratio as where is the basic reproductive ratio for vertical transmission and is the basic reproductive ratio for horizontal transmission. They used the next-generation method and uncertainty analysis to calculate , with a range from 0.037 to 3.743.
A.8. Bluetongue
Gubbins et al. [13] examined bluetongue virus, an insect-borne infectious disease of ruminants in Europe. They used the next-generation method to determine and then Latin Hypercube sampling and partial rank correlation coefficients to assess the effect of parameter variation. Median values for were 0.81 (cattle only), 0.73 (sheep only) and 0.55 (both cattle and sheep), implying that the disease should not persist. However, the corresponding maximum values were 18.77, 15.48, and 12.82, respectively, suggesting that variation in parameters can have an enormous impact on the outcome, including persistence of the disease. The most significant changes in the estimation of were due to the mosquito biting rate, the extrinsic vector incubation period, the vector mortality rate, the probability of transmission from host to vector, and the ratio of vectors to cattle and sheep.
A.9. Plant Pathogens
van den Bosch et al. [6] described a systematic method to calculate the basic reproduction ratio from knowledge of a pathogen's life cycle and its interactions with the host plant. They developed a system of linear difference equations and rearranged the dominant eigenvalue to find . A fraction of infected spores are deposited in location 1, while a fraction are deposited in location 2. This results in where the first term represents the number of spores that successfully germinate in location 1 and the second term represents the number of spores that successfully germinate in location 2. Here, is the probability that a spore deposited on location will germinate, is the number of spores produced per time unit on location , is the infectious period of a spore-producing lesion on location , is the probability that a spore is deposited on a susceptible site, and is the density of susceptible sites in a host population.
They included a separate box explaining pitfalls in calculating the basic reproductive ratio from nonlinear models and noted its nonuniqueness, calculating two distinct values for . However, they were unable to decide which value was correct. The authors conclude by stating that βAlthough nonlinear differential equations are invaluable tools for studying epidemics, they should not serve as the primary basis for the determination of .β
Cunniffe and Gilligan [62] considered the effects of host demography upon the outcomes for invasion, persistence and control of pathogens in a model for botanical epidemics. They used the next-generation method to calculate partitioned into primary and secondary infection. These values are unaffected by host demography. However, the endemic level of infection is highly sensitive to changes in when . If is increased by shortening the infectious period, then the endemic level of infection increases monotonically. However, if increases in transmission rates or decreases in the decay of free-living inoculum increase , then the endemic level of infection first increases (for ) but then decreases (for ). Consequently, increasing the intensity of control can result in more endemic infection.
Acknowledgments
The authors thank Bernhard Konrad, Elissa Schwartz, Jane Heffernan, Lindi Wahl, Helen Kang, and Romulus Breban for technical discussions and are grateful to an anonymous reviewer for comments that significantly improved the paper. R. J. Smith? is supported by an NSERC Discovery Grant, an Early Researcher Award, and funding from MITACS. For citation purposes, note that the question mark in Smith? is part of his name.