Reproduction numbers of infectious disease models.

This primer article focuses on the basic reproduction number, ℛ0 , for infectious diseases, and other reproduction numbers related to ℛ0 that are useful in guiding control strategies. Beginning with a simple population model, the concept is developed for a threshold value of ℛ0 determining whether or not the disease dies out. The next generation matrix method of calculating ℛ0 in a compartmental model is described and illustrated. To address control strategies, type and target reproduction numbers are defined, as well as sensitivity and elasticity indices. These theoretical ideas are then applied to models that are formulated for West Nile virus in birds (a vector-borne disease), cholera in humans (a disease with two transmission pathways), anthrax in animals (a disease that can be spread by dead carcasses and spores), and Zika in humans (spread by mosquitoes and sexual contacts). Some parameter values from literature data are used to illustrate the results. Finally, references for other ways to calculate ℛ0 are given. These are useful for more complicated models that, for example, take account of variations in environmental fluctuation or stochasticity.

Introduction

Infectious diseases continue to debilitate and to cause death in humans and animals, with new disease-causing pathogens emerging and old pathogens reemerging or evolving. For example, viruses give rise to influenza, measles and West Nile virus, bacteria give rise to anthrax, salmonella, chlamydia and cholera, and protozoa give rise to malaria and trypanosomiasis (sleeping sickness). Disease may be passed directly from person to person by respiratory droplets (e.g., measles), via body secretions (e.g., chlamydia), by biting tsetse flies (e.g., trypanosomiasis) or mosquitoes (e.g., malaria), or by ingestion in food or water (e.g., cholera). Some diseases can be controlled by vaccines, antibiotics, antiviral drugs, reduction in vector populations, increased sanitation or behavioral changes. In order to consider control strategies for a particular disease, it is essential to know features of the pathogen, the mode of transmission and other epidemiological details, since as indicated by the above examples, these differ greatly between diseases.

Mathematical modelling can play an important role in helping to quantify possible disease control strategies by focusing on the important aspects of a disease, determining threshold quantities for disease survival, and evaluating the effect of particular control strategies. A very important threshold quantity is the basic reproduction number, sometimes called the basic reproductive number or basic reproductive ratio (Heffernan, Smith, & Wahl, 2005), which is usually denoted by . The epidemiological definition of is the average number of secondary cases produced by one infected individual introduced into a population of susceptible individuals, where an infected individual has acquired the disease, and susceptible individuals are healthy but can acquire the disease.

In reality, the value of for a specific disease depends on many variables, such as location and density of population. For a few specific diseases, Table 1 gives estimates of gleaned from data in the literature.

Table 1

Estimated Mean Values of from Data.

Disease outbreak and location	ℛ_0	Reference
Smallpox in Indian subcont. (1968–73)	4.5	Anderson and May (1991)
Poliomyelitis in Europe (1955–60)	6	Anderson and May (1991)
Measles in Ghana (1960–68)	14.5	Anderson and May (1991)
SARS epidemic in (2002–03)	3.5	Gumel et al. (2004)
1918 Spanish influenza in Geneva
Spring wave	1.5	Chowell, Ammon, Hengartner, and Hyman (2006)
Fall wave	3.8	Chowell et al. (2006)
H2N2 influenza pandemic in US (1957)	1.68	Longini, Halloran, Nizam, and Yang (2004)
H1N1 influenza in South Africa (2009)	1.33	White, Archer, and Pagano (2013)
Ebola in Guinea (2014)	1.51	Althaus (2014)
Zika in South America (2015–16)	2.06	Gao et al. (2016)

The aim of this review is to elaborate on mathematical ways of finding for ODE disease models in a population, bearing in mind the epidemiological meaning of , and to demonstrate how this and other reproduction numbers can be used to guide control strategies. Section 2 introduces simple models that establish notation and serve as background for later sections. The next generation matrix method to theoretically calculate for ODE models is presented in Section 3. In Section 4, the use of and other reproduction numbers to guide control strategies is shown by defining elasticity indices, and type and target reproduction numbers. Sections 5, 6, 7 apply these ideas to models specific for West Nile virus in birds, cholera in humans, and anthrax in animals, respectively. As suggested by a referee, Zika models are briefly discussed in Section 8. For these diseases, numerical values for model parameters are taken from the literature, with references given for these and for proofs (which are not detailed here). A final section, Section 9, gives references to other approaches for calculating , in particular for models formulated in other ways. Inevitably the reference list is incomplete as there have recently been many articles on infectious diseases (it has been said that there is an epidemic of disease models), many of which determine a basic or control reproduction number.

Simple compartmental disease models

SIR epidemic model

To begin with a simple model, assume that each member of a population is either susceptible, infectious (infected with the disease) or recovered from the disease with life-long immunity. If the disease is short lived compared with the population lifetime, then demography can be ignored. Such a model may be appropriate as a very simple model for seasonal influenza, ignoring features such as immunity from past infections. Let denote the number of susceptible, infectious, recovered individuals at time t. Transmission of influenza is airborne or by respiratory secretions on hands, so this is often modeled by mass action, namely a term , where β is the disease transmission rate constant and is the force of infection. Let denote the mean infectious time (about 5 days for influenza), thus is the recovery rate, and let f denote the fraction of infectious individuals who recover from the disease (thus the fraction die from the disease). The flow diagram for the disease dynamics with compartments S, I and R is given in Fig. 1.

Fig. 1

Flow diagram for the model.

Ordinary differential equations (ODEs) for this model are given by

Initially with , and . There is a disease free equilibrium (DFE) with . Focusing on the I equation, the initial behavior is governed by the sign of , or equivalently . This leads to the definition of , with the DFE locally asymptotically stable (LAS) if , but unstable if . This is the product of the transmission rate, the mean infectious time and , and clearly fits with the epidemiological definition of given in the Introduction. Note that is independent of the fraction dying from the disease. From the dynamics of the system, if , then the number of infectious individuals decreases monotonically to 0; whereas if , then this number first increases (before tending to zero); thus acts as a sharp threshold between the disease dying out or causing an epidemic.

If , then initially is approximately . With a knowledge of the mean infectious time, it may be possible to use this formula with initial data from an epidemic to estimate the value of . If and (no death due to disease), then the fraction of susceptibles at the end of the epidemic, where N is the total population, is given by the unique root between 0 and 1 of the final size equation

If data on this fraction are available, then this equation can be used to estimate after the epidemic has passed.

If a vaccine that is perfect is available, and a fraction p of the population is vaccinated, then the disease will not spread if , that is a fraction is vaccinated, giving herd immunity. So this simple model predicts that diseases with large values of require a large proportion of the population to be successfully vaccinated (compare values for smallpox and measles as given in Table 1). Interestingly, the disease that has historically caused the most human deaths, namely smallpox, is the only human disease that currently has been declared eradicated (due largely to vaccination campaigns).

Simple demography can be included in this model with denoting the input of individuals (all susceptible) per unit time, and denoting the natural death rate. The equation for the infectious individuals then becomesgiving with .

Compartmental model

In many infectious diseases there is an exposed period also called a latent period after transmission of infection but before the infected individual can transmit the infection. During this time the pathogen is in the host, but in low numbers so that the host is not yet infectious. If the exposed period is relatively long, then an exposed compartment should be included to give an model. Let E denote the number of exposed individuals, and the mean exposed period be where is the rate of loss of latency. To include simple demography, let A be the input of individuals per unit time, and d be the natural death rate. Assuming that the disease does not cause death, the disease evolves according to the equationswith nonnegative initial conditions. The flow for this model is depicted in Fig. 2. The DFE is given by . There are now two infected compartments (E and I), and linearizing about the DFE it can be seen that only the equations in these two variables determine the stability of the DFE. From these two equations, the Jacobian matrix at the DFE gives the characteristic equation

Fig. 2

Flowchart for the model.

This equation has all roots with negative real parts (thus the DFE is linearly stable) if and only if each coefficient is positive, i.e., . Writing this condition in non-dimensional form as is more epidemiologically meaningful. In fact is the fraction of individuals progressing from exposed to infectious, and is the average infectious time taking death into account. From the biological definition,for this model. Note that if this simplifies to the for the model obtained in Section 2.1. If , then the disease dies out; whereas if , the quadratic has a positive root, and the disease invades the population. In fact the disease then tends to an endemic state with . This endemic equilibrium exists only for , and there is a forward bifurcation at . A mathematical way to determine for this and more complex models is given in the next Section.

Computation of by the next generation matrix method

The Jacobian method used for the model yields a biologically reasonable , but for more complex compartmental models, especially those with more infected compartments, the method is hard to apply as it relies on the algebraic Routh-Hurwitz conditions for stability of the Jacobian matrix. An alternative method proposed by Diekmann, Heesterbeek, and Metz (1990) and elaborated by van den Driessche and Watmough (2002) gives a way of determining for an ODE compartmental model by using the next generation matrix. Here an outline of this method is given, the proofs and further details can be found in van den Driessche and Watmough (2002) and van den Driessche and Watmough (2008).

Let be the number of individuals in each compartment, where the first compartments contain infected individuals. Assume that the DFE exists and is stable in the absence of disease, and that the linearized equations for at the DFE decouple from the other equations. The assumptions are given in more details in the references cited above. Consider these equations written in the form for . In this splitting, is the rate of appearance of new infections in compartment i, and is the rate of other transitions between compartment i and other infected compartments. It is assumed that and , and if .

Now define and for From the biological meanings of F and V, it follows that F is entrywise non-negative and V is a non-singular M-matrix (see Berman and Plemmons (1994)), so is entrywise non-negative. Let be the number of initially infected individuals. Then is an entrywise non-negative vector giving the expected number of new infections. Matrix has entry equal to the expected number of secondary infections in compartment i produced by an infected individual introduced in compartment j. Thus is the next generation matrix andwhere denotes the spectral radius. The linear stability of the DFE is determined from the Jacobian matrix by , where s denotes the maximum real part of the eigenvalues sometimes called the spectral bound. Using the above notation, the relation between this quantity and is given in the following result, the proof of which uses properties of M-matrices; see the references cited above.

If is a DFE of the system , then is locally asymptotically stable if , but unstable if , i.e. sign sign

This next generation matrix approach is now illustrated by returning to the Model of Section 2.2. The infected compartments are E and I. At the DFE matrices F and V aregiving

So has eigenvalues 0 and whereas deduced biologically in Section 2.2. Here is the infection rate in a population of susceptibles, is the fraction progressing from E to I, is the mean time in I, thus the entry of is the expected number of secondary infections produced in compartment E by an infected person originally in E.

As an extension of this model, suppose that individuals in E are mildly infectious at a reduced rate with . This could also be thought of as an initial infectious class before the pathogen has fully developed in the host. Proceeding as previously yieldswith V unchanged. Thus

The first term gives contributions from the infectious compartment I as before, whereas the second term gives contributions from the mildly infectious E compartment, thus increasing the value of .

As another extension of the model (with the E compartment assumed not infectious), suppose that a fraction p of individuals are vaccinated at recruitment into the susceptible population. This is an approximation for vaccination of babies against childhood diseases. In addition, suppose that the vaccine is perfectly effective, so everyone receiving the vaccine is protected from the disease. With input , keeping the total population constant, enter S and enter the R compartment. The DFE then becomes and the vaccination reproduction number also called a control reproduction number, denoted here by , is given by

This gives the fraction that need to be vaccinated as (with ) to bring below the threshold value of one, as determined to give herd immunity (see Section 2.1).

The examples given above each have the rank of the next generation matrix equal to one, so the reproduction number is easy to calculate as the trace of this matrix. A more complicated model in which this occurs is when the next generation matrix is separable, i.e., has entry . Assume that and with as given in Diekmann, Heesterbeek, and Britton (2012, Exercise 7.18) for a sexually transmitted disease in groups, with populations given by and the total population normalized to 1. Then , where denotes the mean value of c and denotes the variance of c.

As a cautionary remark, note that Theorem 1 assures local stability if . For many models, global stability can be proved by comparison theorems (see, for example, Castillo-Chávez, Feng, and Huang (2002)) or Lyapunov functions (see, for example, Shuai and van den Driessche (2013)). However, in some models that have more complicated structure (e.g., include vaccination that is not perfect (Arino et al., 2003, Brauer, 2004, Kribs-Zaleta and Velasco-Hernandez, 2000)) multigroup models (Hadeler & Castillo-Chávez, 1995) endemic equilibria may exist near the DFE for , and there is a backward bifurcation at , with the DFE locally but not globally stable in a range of values below 1. In such a situation, the initial numbers in the infected compartments determine whether or not the disease persists; see for example, Castillo-Chávez and Song (2004) and Martcheva (2015, Section 7.5). So for disease control, may have to be decreased further below 1.

Use of reproduction numbers to suggest control strategies

Related reproduction numbers, and elasticity indices of are now defined, and then in later sections applied to specific disease models.

Type and target reproduction numbers

Herd immunity is applicable if a control strategy is aimed at all individuals in a population, but if the population has several host types, then a control strategy may be aimed at only one host population. For example, in a vector-host model for malaria, spraying may be targeted at the mosquito vectors. To address this Roberts and Heesterbeek (2003) and Heesterbeek and Roberts (2007) introduced the concept of a type reproduction number denoted by . In general such a strategy influences one row or one column of the next generation matrix, depending on whether the control influences susceptibility or infectivity. This concept was further generalized by Shuai, Heesterbeek, and van den Driessche (2013) by singling out sets of entries for control, leading to the definition of a target reproduction number. For example, in a cholera outbreak, the prevention of contacts among children could be a reduction strategy; an age structured model would then target a specific entry of the next generation matrix.

To describe this algebraically, let be the order n next generation matrix (i.e., ) and suppose that the entries in a set S are to be targeted with a control strategy (this can be either a decrease or an increase in entries of S). The target matrix has entry equal to if entry , and 0 otherwise. If , then the target reproduction number is defined aswhere is the identity matrix of order n. When all entries in one or more rows (or columns) of K are targeted, then the target reproduction number reduces to the type reproduction number as defined by Roberts and Heesterbeek (2003) and Heesterbeek and Roberts (2007). If a fraction at least of the contributions due to the set S can be prevented, then the disease will die out. Note that if K is irreducible, then if and only if (Roberts and Heesterbeek, 2003, Shuai et al., 2013).

Sensitivity and elasticity

To determine best control measures, knowledge of the relative importance of the different factors responsible for transmission is useful. Initially disease transmission is related to and sensitivity predicts which parameters have a high impact on . The sensitivity index of with respect to a parameter ω is Another measure is the elasticity index (normalized sensitivity index) that measures the relative change of with respect to ω, denoted by , and defined as

The sign of the elasticity index tells whether increases (positive sign) or decreases (negative sign) with the parameter; whereas the magnitude determines the relative importance of the parameter. These indices can guide control by indicating the most important parameters to target, although feasibility and cost play a role in practical control strategy. If is known explicitly, then the elasticity index for each parameter can be computed explicitly, and evaluated for a given set of parameters. The magnitude of the elasticity indices depends on these parameter values, which are probably only estimates. To help identify the robustness of to the parameters, Latin Hypercube Sampling maximin criteria can be used; see, for example, Blower and Dowlatabadi (1994). Another technique to investigate this is to compute over the feasible region of a given parameter while keeping the other parameters fixed at baseline values; see for example Manore, Hickmann, Xu, and Hyman (2014).

A model for West Nile virus

Some diseases, for example, West Nile virus, dengue fever, malaria, Zika virus are transmitted through a vector (for these diseases the vectors are various species of mosquitoes), rather than directly from person to person. Female mosquitoes bite to obtain a blood meal that is essential for reproduction, so only female mosquitoes need be considered. West Nile virus can kill birds and humans, but infected mosquitoes remain infectious for life and do not die of the virus. Birds can transmit the virus back to mosquitoes, whereas humans appear to be dead end hosts, and so may be excluded from a simple model. Since the life-cycle of mosquitoes is much shorter than that of birds, mosquito demography should be included in a model, but bird demography can be ignored. For a simple mosquito-bird (vector-host) model, take S and I compartments in each population, giving a system of four ODEs as formulated by Wonham and Lewis (2008). The equations for this system with nonnegative initial conditions arewhere the variables and parameters are defined as

: number of susceptible birds, mosquitoes

: number of infectious birds, mosquitoes

: probability per bite of virus transmission to bird, mosquito

β: biting rate of mosquitoes on birds

: bird death rate from virus

: mosquito birth, natural death rate

The disease is transferred by an infectious mosquito biting a susceptible bird, or by a susceptible mosquito biting an infectious bird, This cross infection between mosquitoes and birds is illustrated in the flow diagram in Fig. 3. The transmission is assumed to be frequency dependent; see Wonham and Lewis (2008) for more discussion on the model and transmission assumptions.

Fig. 3

Flowchart for the West Nile virus model by Wonham and Lewis (2008).

Assuming that so the bird population is constant, there is a DFE with all birds and mosquitoes susceptible, with populations denoted by , and . Using the next generation matrix at the DFEgiving

The first ratio under the square root represents the number of bird infections caused by one infectious mosquito, and the second represents the number of mosquito infections caused by one infected bird. The square root represents a geometric mean. In the literature the square root is often omitted, giving the same threshold for stability at 1, but taking the average number of secondary infected humans resulting from a single infected human; see, for example, Roberts and Heesterbeek (2003).

The relative importance of each parameter for control can be estimated by computing elasticity indices. For this model, . Thus reducing the biting rate of mosquitoes has the largest proportional effect on reducing .

Targeting the mosquito to bird transmission, the target reproduction number has , thus

If the transmission from mosquitoes to birds can be reduced by a fraction at least , then the vector-host disease will die out.

This simple model has in particular neglected the period during which mosquitoes are in the larval stage, and also the exposed period of infected mosquitoes during which the viral load becomes sufficiently high for bites to be able to transmit the disease. Both these periods are significant fractions of the mosquito life-span. Including the two extra compartments of larval and exposed mosquitoes, gives a 6-dimensional ODE system, as formulated in Wonham and Lewis (2008). The resulting is not changed by the larval stage, but the exposed class reduces by a factor that is the square root of the probability of an exposed mosquito becoming infectious. Numerical model simulations in Wonham and Lewis (2008) show that the inclusion of these classes delays the disease outbreak.

Jiang, Qiu, Wu, and Zhu (2009) consider a bird-mosquito West Nile virus model (proposed by Bowman, Gumel, van den Driessche, Wu, and Zhu (2005)) and find that backward bifurcation is possible for some parameter values, thus the initial numbers of mosquitoes and birds are important in determining whether or not the disease dies out, even if . More complicated vector-host models require appropriate parameter values to estimate elasticity indices and so to guide control planning; see, for example, Manore et al. (2014) for a comparison of dengue and chikungunya dynamics, and Cai, Li, Tuncer, Martcheva, and Lashari (2017) for a malaria model that may also exhibit backward bifurcation. Temperature plays an important role in the reproduction of many vectors that transmit vector-borne diseases. Thus climate change and variations in temperature may have an important influence on and disease persistence. A detailed study of this effect for Lyme Disease, which is transmitted by ticks, using data from northeastern North America is given by Ogden et al. (2014), and concludes that climate warming may have partly driven the emergence of Lyme disease in this region. More recently, Wang and Zhao (2017) develop a periodic time-delayed model of Lyme disease, compute from data, and show that can be driven below 1 if the recruitment rate of tick larvae is reduced.

Cholera models

Cholera is an infection of the small intestine caused by the bacterium Vibrio cholerae. Infection causes mild diarrhea in most cases, but some cases develop severe diarrhea and vomiting, which if untreated may lead to death within a few hours due to dehydration and electrolyte imbalance. Cholera can be transmitted indirectly to humans via water infected with the bacterium that is shed from infected humans, or directly from human to human (without differentiating between males and females). The relative importance of these two transmission paths is a key factor in designing control strategies as illustrated by the following models of cholera from the literature.

Direct and indirect transmission, Tien and Earn (2010)

The same notation is used as previously, augmented by a variable B measuring the concentration of bacterium in the water, with δ the removal rate of the bacterium, and ξ the shedding rate of bacterium by infectious humans. The model formulated by Tien and Earn (2010) incorporates both direct and indirect transmission by using mass action terms with transmission rate constants and . Their model assumes that there is no death due to cholera, so that the human population remains constant. The flow for this model of cholera dynamics is given in Fig. 4. The equations governing the dynamics arewith nonnegative initial conditions having and . The DFE is given by , and assuming that shedding is not a new infection (i.e., goes into the V matrix), the next generation matrix has rank 1 and gives

Fig. 4

Flowchart of the cholera model by Tien and Earn (2010).

The two routes of transmission enter in an additive way, thus , with and . If either the indirect or direct reproduction number is greater than one, then , emphasising the fact that both transmission routes must be controlled to eliminate cholera.

However, if shedding is regarded as a new infection, then this alternative splitting gives

This splitting gives a next generation matrix of rank 2, so is the positive root of the quadratic equation

From the signs of the coefficients, this quadratic has a unique positive root, ; and if , then . But , giving the same threshold as derived for above. However, it is hard to interpret biologically.

For qualitative control, both splittings suggest vaccination to reduce the susceptible population at the DFE, however current vaccines are only effective for a short duration. More importantly, provision of clean water, which reduces indirect transmission (decreases ), water treatment to increase the bacterium decay rate δ, improvement in sanitation by disposal of human faeces (decreases ξ), and improvement in hygiene (decreases ), all combined in the acronym WASH, may be more beneficial in controlling cholera.

For this cholera model with shedding not a new infection, elasticity indices of can be calculated explicitly as

The parameter estimates for the 2006 cholera outbreak in Angola as reported in Table 1 of Eisenberg, Robertson, and Tien (2013) with time unit one day, namely , give

This latter result means that if the recovery rate is increased, then decreases by the same relative amount. From these elasticity indices for the normalized transmission rates, it follows that the indirect transmission has the larger effect on . Note that depends on the water transmission rate as well as the shedding rate and pathogen decay rate. However, Eisenberg et al. (2013) find that both transmission pathways are needed to explain the data well, and both should be considered in formulating public health policies to reduce the impact of cholera.

Taking the next generation matrix , type and target reproduction numbers can be calculated; see Shuai et al. (2013). If direct transmission is targeted (for example, by improving hygiene or isolation), then , and the formula in Section 4.1 gives this target reproduction number as provided that , i.e., cholera cannot be maintained by only indirect transmission. If adequate clean water is provided, then the target reproduction number with gives provided that cholera cannot be maintained by direct transmission alone. Thus if clean water is provided to at least a proportion of the population, then cholera will die out.

A random network model, Li, Ma, and van den Driessche (2015).

To take account of realistic heterogeneity in human contacts, network models have been introduced; see, for example, Newman (2003). The network consists of nodes as humans and edges as contacts. Miller (2011) showed that an model on a configuration type random network can be represented by a 2-dimensional ODE system. Li et al. (2015) extended this model to a random network of humans with each connected to a single infectious water node, as is appropriate for cholera.

Assume that the human to human network is generated by the probability generating function where is the degree distribution and k is the node degree. Let be the transmission rate along a human to human edge, and let with the rate of transmission along the water to human edges caused by a unit concentration of bacterium in the water. Parameters are as in the previous model (Section 6.1). The dynamics of this network cholera model are then given by a 3-dimensional ODE system, see Li et al. (2015) for details, which has the next generation matrix at the DFE given bywhere , and are the basic reproduction numbers for the human to water to human, the water infected human to human, and the human infected human to human infected individuals, respectively. Note that is the average degree of the network. Since has rank 2 (the second row and column can be ignored in calculating the spectral radius) the basic reproduction number is given by the largest positive root of the quadratic equationwith the disease threshold given by . For a Poisson network, the variance is equal to the mean, i.e., , giving , thus as in the compartmental model of Section 6.1. However, for other networks this additive result of indirect and direct transmission does not hold. For example, on a scale-free network, the variance is greater than the mean, so and . Thus the network structure is important in considering ways to control cholera by bringing below the threshold value of 1.

If , then target reproduction numbers can help guide strategies for control. For example, if , i.e., cholera cannot be sustained in the population just by the human to human transmission, then the type reproduction on the water transmission is given by

Preventing a fraction of at least contributions from human infections due to water infection would cause cholera to die out.

For a Poisson network, the elasticity indices with respect to the transmission parameters are

For other networks, the elasticity indices can be found by implicitly differentiating the quadratic in z.

A model for anthrax transmission

Anthrax is caused by a spore forming bacterium Bacillus anthracis. It is mainly a disease of animals, but humans can become infected by contact with infected animals, their hides, wool or meat. Infected animals may recover from anthrax (usually carnivores) or die from the disease (usually herbivores) (WHO, 2016). Transmission to animals may occur from eating infected carcasses or by inhalation of spores, and rarely by direct transmission. A model including these first two important models of transmission in animals can be formulated as a system of four ODEs, and is a special case of the model in Saad-Roy, van den Driessche, and Yakubu (2017) where direct transmission is also included. Let denote the number of susceptible animals, infected animals and infected carcasses, and A denote the grams of anthrax spores in the environment. It is assumed that animals follow logistic growth with birth rate r and carrying capacity K, with all newborn animals being susceptible. Spores grow at rate β per carcass and decay at rate α. Transmission rates from infected carcasses and from spores are denoted by and , respectively. Parameters and κ denote the animal natural death rate, death rate due to anthrax, recovery rate of infected animals, carcass consumption rate and carcass decay rate, respectively. The dynamics of the model are given by the system

Assuming that , the total animal population persists and there is a DFE with . Taking the infected compartments as , and C and assuming that the transmission terms, spore growth rate and death of infected animals give rise to new infections (as in Friedman and Yakubu (2013)), the Jacobian matrix at the DFE is decomposed as where is a nonnegative matrix and is a nonsingular M-matrix. Thus,

This gives as the positive root of the cubic polynomialwhere and .

Alternatively, taking spore growth rate and death of infected animals as transfer, the Jacobian is split in a different way, givingand , which is the sum of infections from spores in the environment () and from feeding on carcasses . Note that , showing that if and only if . Thus gives the same threshold as the more biologically meaningful basic reproduction number for anthrax to die out or persist.

For values of the growth parameters and giving , Table 2 gives baseline parameter values from Saad-Roy et al. (2017), elasticity indices of and numerical values for these indices. These baseline parameters give , , thus . Note that the model in Saad-Roy et al. (2017) includes a small amount of direct transmission between animals. Also in the caption of Table 2 in Saad-Roy et al. (2017) the values of and are erroneously switched: the correct values given here result in switching of the elasticity indices of and , and the elasticity index of α should be in Saad-Roy et al. (2017).

Table 2

Baseline Parameter Values and Elasticity Indices for the Anthrax Model.

Parameter	Value	Elasticity	Numerical Elasticity
δ	0.05	ϒδℛ_0=−δS_0δS_0+κ	−0.962
η_c	0.1	ϒη_cℛ_0=ℛ_0cℛ_0	0.909
τ	0.1	ϒτℛ_0=−τγ+d+τ	−0.409
γ	0.143	ϒγℛ_0=−γτ(γ+d+τ)(γ+d)	0.404
η_a	0.5	ϒη_aℛ_0=ℛ_0aℛ_0	0.091
β	0.5	ϒη_aℛ_0=ℛ_0aℛ_0	0.091
α	0.1	ϒαℛ_0=−ℛ_0aℛ_0	−0.091
κ	0.1	ϒκℛ_0=−κδS_0+κ	−0.038

For these baseline parameter values, the elasticity indices show that the carcass decay rate and transmission rate from infected carcasses have the largest relative influence on . Thus a control strategy that increases carcass decay rate or decreases carcass transmission is likely to be successful. This can be further explored by calculating the type reproduction number focussing on the removal of infected carcasses on death. Taking in , gives . Thus for eradication at least a fraction of infected carcasses need to be removed, which is about for the baseline parameters in Table 2. Further discussion on removal of carcasses at any time after death is given in Sections 5.1 and 6.1 of Saad-Roy et al. (2017).

Modelling Zika transmission

An outbreak of Zika virus started in Brazil in early 2015, and since then has spread to many parts of the Americas. Although Zika symptoms are usually mild, it is of great concern as serious brain anomalies, such as microcephaly, can occur in fetuses of infected pregnant females. Zika is transmitted to humans by mosquitoes and also sexually by an infected male to a susceptible female. These two modes of transmission are included in models in the recent literature; see, for example, Brauer et al., 2016, Gao et al., 2016, Saad-Roy et al., 2016. Neglecting the exposed compartments, a simplified version of the model in Brauer et al. (2016) is

For this system, the human population is divided into susceptible humans, infectious humans and recovered humans (this equation uncouples); the constant mosquito population is divided into susceptible mosquitoes and infectious mosquitoes. The positive parameters are γ the rate of human recovery, the birth and death rate of mosquitoes, the biting rate times the probability of transmission from mosquito to human, the biting rate times the probability of transmission from human to mosquito, and α the direct sexual transmission from human to human.

Using the next generation matrix method with all infection terms in matrix F, the basic reproduction number is the positive root of the quadratic equation , with and , accounting for the direct sexual transmission and vector transmission, respectively. Note that in this model the transmission term from mosquito to human is normalized differently than that from mosquito to bird in the West Nile virus model of Section 5, and also note similarity between this model and the cholera model in Section 6.1 with direct and indirect transmission. Brauer et al. (2016) show different splittings for matrices F and V leading to different expressions for the basic reproduction number, one such splitting giving . This point about multiple expressions for the basic reproduction number (in both discrete and continuous time models) is elaborated by Cushing and Diekmann (2016).

In their Zika model, Gao et al. (2016) include exposed compartments and also divide the infected human compartment into symptomatically, convalescent and asymptomatically infected. Their global sensitivity analysis shows that their basic reproduction number is most sensitive to mosquito biting and death rates. Saad-Roy et al. (2017) divide the human population in two regions (a source and import region) into sexually actives males and females and sexually inactive humans. Sexually active males are further divided into two stages, since Zika can be transmitted from males to females by blood and by semen, with the virus staying longer in the semen. Sexual contacts between males and females are represented by a random bipartite directed (male to female) network, with homogeneous mixing between mosquitoes and humans. Input parameter values are used to compute from a matrix of order 14, and a value of is found. This model predicts that by the time one microcephaly case is detected, a sizable fraction of the population has been infected with Zika, emphasizing the need for surveillance of this disease.

Other approaches to obtain

This concluding section gives a glimpse of other results about the basic reproduction number, including calculating it for models given by systems other than ODEs, and estimating its magnitude from data. A few relevant references are given so that a reader can follow up on these glimpses to obtain more details.

Next generation operator to calculate

The method of computing as the spectral radius of the next generation matrix given in Section 3 assumes that the population is divided into discrete compartments or states. However, if these are continuous, as for example in an age structured model, then the next generation operator K specifies how the disease evolves in a generation, namely , where is the kernel of the linear positive operator K, and is a subset of ; see Diekmann et al. (1990) and Diekmann et al. (2012, Section 7.3). The kernel gives the expected number of new cases with state x caused by an infectious individual in state y during its period of infectiousness. Then the basic reproduction number is defined as the spectral radius of the positive operator K. If the range of K is one dimensional, then is the unique non-zero eigenvalue; see Diekmann et al. (2012, Section 7.4) for some examples, and see Section 4 in Thieme (2009) for discussion of an model with variable susceptibility.

in a periodic environment

The method of Section 3 assumes that the environment is constant, but in reality environmental parameters are changing, for example temperature changes periodically. The models then become non-autonomous dynamical systems. In the past twenty or so years, many authors have extended the definition of to periodic environments; see for example, Bacaër and Ait Dads, 2011, Inaba, 2012, Thieme, 2009, Wang and Zhao, 2017, Zhao, 2013, Chapter 11) and references therein. Inaba (2012) introduces a new definition of in a heterogeneous environment based on the generation evolution operator as a generalization of previous definitions.

Survival function method to compute

Let denote the infection survival probability, that is the probability that a newly infected individual remains infected (and infectious) for at least time a, and denote the average number of newly infected individuals that an infectious individual produces per unit time. Using these functions, . For the model of Section 2.1, the infectious period is assumed to have a negative exponential distribution. Thus and . Then , as determined in Section 2.1.

This method can be used for models that take account of the fact that infectivity (the probability of transmission from a contact) varies with time since infection, rather than being constant throughout the infectious period. In such an age of infection SIR model, with a the time since first infection, where is the DFE equilibrium. This method is not restricted to models described by ODEs and can be extended to a series of states, for example to a vector-host system; see, for example, Heffernan et al. (2005, Section 2.1).

Calculation of for discrete time systems

There are many discrete time epidemic models in the literature; see, for example, Allen (1994) and Yakubu (2010). Care must be taken in formulating these discrete time models so that the number of individuals in each compartment does not go negative.

The next generation matrix approach can also be used for discrete time epidemic models; see Allen and van den Driessche (2008) and references therein. To describe this briefly, assume that the equations of the epidemic model are given by where denotes the population states at time t, with infected and the remaining states uninfected. Assuming that there is a unique DFE, then linearizing the infected state equations about this DFE gives where F and T are nonnegative matrices evaluated at the DFE. Here F is the matrix of new infections that survive the time interval, and T is the transition matrix with . Then matrix is the next generation matrix for this discrete system, and the basic reproduction number . Allen and van den Driessche (2008) give details of the assumptions and show that the DFE is locally asymptotically stable if but unstable if . The assumed order of events within each time step is important, since different assumptions can yield different values; see Lewis, Rencławowicz, van den Driessche, and Wonham (2006) for an illustration of this on discrete time West Nile virus models.

A graph theoretic method to determine the net reproductive rate for discrete time systems was developed by de-Camino-Beck and Lewis (2008) and applied to the control of invasive species. They use graph reduction rules on the life-cycle graph to calculate . An analogous method for continuous time models given by ODEs is described in de-Camino-Beck, Lewis and van den Driessche (2009), where an algorithm to compute by reduction on the digraph associated with the matrix is given.

Stochastic models

An excellent introduction to stochastic epidemic models formulated as discrete time Markov chains, continuous time Markov chains, and stochastic differential equations is given in Allen (2010, Chapter 3), where some numerical examples of stochastic sample paths and corresponding deterministic solutions are given. Also, Keeling and Rohani (2008, Chapter 6) give three methods for approximating stochasticity into disease transmission and recovery. The probability of an outbreak depends on the value of and also on the initial number of infected individuals, . If the population size is large, then this probability is approximately 0 if , but if . This estimate applies to the stochastic and models only for a finite time range, since as the probability of an outbreak is zero (an absorbing state). Thus, one difference between stochastic and deterministic models is that stochastic sample paths can converge to the disease free state, whereas the corresponding deterministic solution converges to an endemic equilibrium. Thus in a vaccination program, stochasticity may lead to disease extinction before the herd immunity level of vaccination is reached.

Estimation of from data

As described in previous sections, the magnitude of is important in determining the severity of a disease and for designing control strategies. Early in an outbreak of some diseases, for example seasonal influenza, it may be fairly easy to estimate the average length of infection and the initial number of infectious individuals, thus an estimate of can be obtained from the formula given in Section 2.1. If the exponential growth rate of the initial phase of the outbreak, usually denoted by r, is observed from data, then . However, this formula is based on a simple model, and so has limited applicability. A recent paper by Sanches and Massad (2016) compares three methods to estimate for dengue fever, without this assumption of initial exponential growth. Ma, Dushoff, Bolker, and Earn (2014) use maximum likelihood to compare four commonly used phenomenological models (exponential, Richards, logistic, delayed logistic) in estimating initial growth rates of simulated epidemic models. Moreover, in a series of recent papers, Chowell and collaborators assume sub-exponential initial growth given bysee Chowell and Viboud, 2016, Chowell et al., 2016a, Chowell et al., 2016b, Chowell et al., 2016c, and Viboud, Simonsen, and Chowell (2016). They find that the suboptimal growth model outperforms the exponential model (with ) for several infectious disease datasets, and that the effective reproduction number declines within 3–5 disease generations (rather than maintaining a constant ).

Most other methods to estimate require data over a longer time period or until the outbreak is over. For example, after an epidemic when the proportion of susceptibles who did not catch the disease is known (assuming that all individuals are originally susceptible), may be estimated from the final size equation, giving ; see Section 2.1. Using standard optimization techniques, parameters in a model can be estimated from data points at successive times during an outbreak. However, this means that some time must have elapsed for sufficient data to be collected, and so gives a delay in estimating , especially for an emerging disease. In addition, stochastic effects are important, particularly near the beginning of an outbreak. For more details on estimating and fitting models to data, see Chowell and Brauer (2009, Sections 8, 9), Diekmann et al. (2012, Chapter 13), and Martcheva (2015, Chapter 6).