Bulletin of Mathematical Biology - Semantic Scholar

Bulletin of Mathematical Biology (2006) 68: 491–509 DOI 10.1007/s11538-005-9039-7

ORIGINAL ARTICLE

A Comparison of Continuous and Discrete-Time West Nile Virus Models Mark A. Lewisa,b,∗ , Joanna Rencławowicza,c,d , P. van den Driessched , Marjorie Wonhamb a

Department of Mathematical and Statistical Sciences, University of Alberta, Edmonton, Canada AB T6G 2G1 b Department of Biological Sciences, University of Alberta, Edmonton, Canada AB T6G 2G1 c ´ Institute of Mathematics, Polish Academy of Sciences, Sniadeckich 8, 00-956 Warsaw, Poland d Department of Mathematics and Statistics, University of Victoria, Victoria, Canada BC V8W 3P4 Received: 3 February 2005 / Accepted: 8 April 2005 / Published online: 5 April 2006 C Society for Mathematical Biology 2006

Abstract The first recorded North American epidemic of West Nile virus was detected in New York state in 1999, and since then the virus has spread and become established in much of North America. Mathematical models for this vector-transmitted disease with cross-infection between mosquitoes and birds have recently been formulated with the aim of predicting disease dynamics and evaluating possible control methods. We consider discrete and continuous time versions of the West Nile virus models proposed by Wonham et al. [Proc. R. Soc. Lond. B 271:501–507, 2004] and by Thomas and Urena [Math. Comput. Modell. 34:771–781, 2001], and evaluate the basic reproduction number as the spectral radius of the next-generation matrix in each case. The assumptions on mosquitofeeding efficiency are crucial for the basic reproduction number calculation. Differing assumptions lead to the conclusion from one model [Wonham, M.J. et al., Proc. R. Soc. Lond. B 271:501–507, 2004] that a reduction in bird density would exacerbate the epidemic, while the other model [Thomas, D.M., Urena, B., Math. Comput. Modell. 34:771–781, 2001] predicts the opposite: a reduction in bird density would help control the epidemic. Keywords West Nile virus · Basic reproduction number · Discrete time model · Disease control · Next generation operator · Spectral radius

∗ Corresponding

author. E-mail addresses: [email protected] (Mark A. Lewis), [email protected] (Joanna Rencławowicz).

492

Bulletin of Mathematical Biology (2006) 68: 491–509

1. Introduction Although West Nile (WN) virus is endemic in Africa, the Middle East and western Asia, the first recorded North American epidemic of WN virus was detected in New York state in 1999. Since 1999 WN virus has spread and has become established in much of North America. Recent mathematical models for this disease have been proposed in an attempt to predict disease dynamics and elucidate control methods (Bowman et al., 2005; Lord and Day, 2001; Thomas and Urena, 2001; Wonham et al., 2004). The temporal spread of WN virus involves an interplay of the transmission between birds via female mosquitoes as disease vectors and of the disease dynamics within a reservoir of birds. While birds can die quickly from the virus (especially corvids, such as crows and jays), the mosquito disease vectors do not appear to be affected adversely by the disease. The interaction of WN virus with secondary hosts (mainly humans and horses) is dead-end in the sense that there is no evidence that these secondary hosts can infect feeding mosquitoes (although see Higgs et al. (2005) for very recent evidence that some mammals may not be dead-end hosts). Thus a cross-infection model involving susceptible-infectious interactions between birds (reservoirs) and mosquitoes (vectors) is the central starting point for formulating WN virus dynamics, and is the basis for all models considered here. Given a mathematical model for disease spread, the basic reproduction number, R0 , is an essential summary parameter. It is defined as the expected number of secondary cases caused by a single infected individual introduced into an otherwise susceptible population. If R0 is greater than one, a local disease outbreak is possible. Control methods can be designed to bring the control basic reproduction number, which we also denote by R0 , to a value less than one. The dependence of R0 on model parameters can be used (through the modification of parameter values) to evaluate efficacy of such control methods. However, different models for a given disease may not deliver the same R0 . Models for the same disease can be based on differing assumptions about disease transmission and dynamics. These, in turn, result in different R0 expressions. In evaluating the effect of control methods by using R0 , the conclusions that are drawn depend crucially on the model assumptions. Thus, if R0 is to be used to evaluate the effect of control measures on disease outbreak, it is necessary to link the R0 value directly to model assumptions. This is the primary purpose of the paper. Indeed, we will show that the assumptions on mosquito-feeding efficiency are crucial for the basic reproduction number calculation. Differing assumptions lead one model (Wonham et al., 2004) to conclude that a reduction in bird density would exacerbate the epidemic, while the other model (Thomas and Urena, 2001) predicts the opposite: a reduction in bird density would help control the epidemic. Our approach is to consider how both biological and temporal model structure influence R0 . We proceed by example, analyzing, in depth, continuous and discrete-time versions of models by Wonham et al. (2004) and Thomas and Urena (2001). Both consider the North American WN virus epidemic, but are formulated with different biological assumptions and time structures. The first model (Wonham et al., 2004) treats time as a continuous variable, and is an eight-dimensional system of ordinary differential equations for vector and reservoir compartments. For comparison, we also formulate two different discrete-time


493

versions of this model. The second model (Thomas and Urena, 2001) is in discrete time and consists of a nine-dimensional difference equation system for vector, reservoir and human compartments. We also formulate a continuous time version of this model. For the resulting five models, we calculate and compare R0 expressions and consider their differing implications for disease control. The secondary purpose of the paper is to demonstrate, via our sample calculations, how the basic reproduction number can be calculated in a straightforward way for both discrete and continuous time models, even when the models become complicated, such as in the eight- and nine-dimensional cross-infection models for WN virus (Thomas and Urena, 2001; Wonham et al., 2004). 2. Methods to calculate R0 2.1. Continuous time models For an epidemic model with disease compartments that is formulated as an ordinary differential equation system, a precise mathematical definition of R0 is the spectral radius of the next-generation matrix (Diekmann and Heesterbeek, 1999; van den Driessche and Watmough, 2002). Provided that a disease-free equilibrium (DFE) exists, and some other biologically realistic conditions are satisfied, then the next-generation matrix can be determined from the system as follows. Write the equations for the infected compartments only as dx = (F − V)(x), dt where vector x gives the number in each infected compartment, F denotes the rate of new infections and V denotes the rate of transfer (by other means) between compartments. Let F and V denote the linearized matrices at the DFE from F and V, respectively. Then F V −1 is the next-generation matrix, with the (i, j) entry giving the expected number of new infections in compartment i produced by an infected individual introduced into compartment j. The basic reproduction number R0 is the spectral radius of F V −1 (see Diekmann and Heesterbeek, 1999; van den Driessche and Watmough, 2002). The DFE is locally asymptotically stable if the matrix F − V has all eigenvalues with negative real parts. With the above definition of R0 , this stability condition can be shown by using M-matrix theory to be equivalent to R0 < 1. In addition, F − V is unstable if R0 > 1. Therefore, if introduced at a low level, the disease dies out when R0 < 1, but persists in the population when R0 > 1. The exact form of R0 is thus important in determining control strategies for the disease. 2.2. Discrete-time models For discrete-time epidemic models, the equations for the infected compartments are written as x(t + 1) = (F + T )x(t) where x(t) is the number in each infected compartment at time step t, F represents the new infections and T represents other transitions between compartments.

494


Linearization at the DFE (which is assumed to exist), gives rise to nonnegative matrices F and T, with the spectral radius of T less than 1. As in Caswell (2001) and Li and Schneider (2002), the discrete next generation matrix, which projects the infected compartments from one-time step to the next, is given by F(I − T)−1 , where I denotes an identity matrix. For a discrete system, R0 is the spectral radius of F(I − T)−1 . It follows from Perron–Frobenius theory (Caswell, 2001; Li and Schneider, 2002), that the DFE is linearly stable or unstable according to whether R0 is less than or greater than one. Thus the exact formulation of a discrete model, which in turn gives an expression for R0 , is important in determining whether or not the disease can persist, and in suggesting control measures. 2.3. Common notation To compare different models, we introduce common notation for all state variables (numbers of mosquitoes, birds and humans in the different disease compartments) and parameters in the two original models. Exposed and infective compartments must both be considered as infected when calculating R0 . Since West Nile virus is a vector-transmitted disease, we expect a square root in the expression for R0 , which arises as a geometric mean of the vector and reservoir variables (Diekmann and Heesterbeek, 1999; van den Driessche and Watmough, 2002). Vector State variables Larval Susceptible Exposed Infectious Recovered Dead Total adults Total Parameters Birth Proportion of births that are infected Maturation Death (natural) Death (from virus) Vector biting on host Virus transmission (to) Virus incubation Recovery from virus

LV SV EV IV

AV NV bV ρV mV dL , dV αV κV

Reservoir

Humans

SR

SH

IR RR XR

IH RH

NR

NH

bR

bH

dR δR βR αR

βH

γR

γH

For continuous time models, all parameters are per capita rates per unit time, except for the proportion ρV and the probabilities αV and αR . For discrete-time models, all parameters except ρV are probabilities or numbers per unit time. In the


495

analysis that follows a superscript ˆ is used to distinguish parameters for discretetime models from parameters for continuous time models. The models we consider here are for a single season under constant environmental conditions (but see Wonham et al., 2004, for extensions to a variable environment). Hence there is birth and death of mosquitoes (vectors), but only death of birds (reservoirs), as bird reproduction is assumed to have occurred before the mosquito season. The Wonham et al. (2004) model allows for virus-induced death of birds (which is necessary for birds, such as corvids, that die quickly from WN virus), while the Thomas and Urena (2001) model does not allow for virus-induced death of birds (which may be a good approximation for birds, such as passerines, that do not die quickly from WN virus) but includes natural death in all bird compartments. The Thomas and Urena (2001) model also allows for vertical transmission of the virus from mosquito parent to offspring, while the Wonham et al. (2004) model does not. The two models differ in their assumptions on the disease transmission terms. The Wonham et al. (2004) model was formulated in continuous time whereas the Thomas and Urena (2001) model was formulated in discrete time, with time steps of a week. The calculation for the basic reproduction number was made for the Wonham et al. (2004) but not for the Thomas and Urena (2001) model. We formulate discrete- and continuous-time versions for both models, and make the calculation of basic reproduction number for each formulation, so as to facilitate model comparison. 3. West Nile virus model of Wonham et al. (2004) 3.1. Continuous time model We now give the continuous time WN virus model as formulated by Wonham et al. (2004). The authors include age structure for the female mosquito population by dividing this population into larvae and adults, with birth into the larval stage, and natural death in each stage. The adult stage is divided into susceptible, exposed (latent) and infectious compartments. For the one season model, they assume that WN virus can cause reservoir death, but natural birth and death of the reservoir population is ignored. This population is divided into susceptible, infectious, recovered and dead compartments. Cross-infection between mosquitoes and birds is modeled by mass action incidence normalized by the total population of birds. This arises since female mosquitoes only take a fixed number of blood meals per unit time, and follows a similar term used to model malaria; see, for example, Anderson and May (1991). Since humans are dead-end hosts, they are not included in this model. In the common notation, the dynamics are given by the following ordinary differential equation system: Vectors (V): dLV = bV (SV + EV + IV ) − mV LV − dL LV dt dSV IR = −αV βR SV + mV LV − dV SV dt NR

496


dEV IR SV − (κV + dV )EV = αV βR dt NR dIV = κV EV − dV IV . dt Reservoirs (R): dSR SR IV = −αR βR dt NR dIR SR IV − (δR + γR )IR = αR βR dt NR dRR = γR IR dt dXR = δR IR . dt For the existence of a disease-free equilibrium it is assumed that vector birth and death rates balance in the absence of disease. This is expressed by the following parameter constraint: bV = dV (1 + dL /mV ). The disease-free equilibrium is (LV , SV , EV , IV , SR , IR , RR , XR ) =

bV A∗V , A∗V , 0, 0, NR∗ , 0, 0, 0 . mV + dL

The infected variables are (EV , IV , IR ) and with F, the rate of appearance of new infections, and V, the rate of transfer between compartments, 

 IR   S V (κV + dV )EV EV   N R                IV  = F − V =  0   −  dV (IV − κV )EV   .            S R IR t (δR + γR )IR αR βR IV NR 



αV βR

The corresponding linearized matrices at the DFE are 

A∗V  0 0 αV βR N∗ R   F = 0 0 0   0 0 αR βR

    ,   



κV + dV 0

0

dV

0

 V =  −κV 0

0 δR + γR

  .


497

Then 

F V −1

    =    

0

0

0

0

αR βR κV αR βR (κV + dV )dV dV

 αV βR A∗V (δR + γR )NR∗      0     0

so the spectral radius of F V −1 is R0 =

A∗V αV αR βR2 κV dV (dV + κV )(δR + γR ) NR∗

as found in Wonham et al. (2004). It can be seen that R0 is the geometric mean A∗V αR βR V βR κV . The first term is the product of the infection rate to of dVα(d ∗ and (δR +γR ) V +κV ) NR the vector at the DFE, the average time that a vector spends in the infective class and the probability that a vector entering the exposed class survives to become infective. The second term is the product of the infection rate to the reservoir and the average time that a bird spends in the infective class before dying or recovering.

3.2. Discrete-time version of Wonham et al. (2004) model—version 1 To uniquely write down the difference equations from the original continuous time model, the ordering of events needs to be specified. Thus, we consider the following set of assumptions. 1. Birth, infection and transfer between compartments occur at the beginning of the time step. 2. Natural and disease-induced mortality occur at the end of the time step. Note that a different ordering of events is considered in the Appendix. Under assumptions (1) and (2), the difference equation system takes the form: LV (t + 1) = (1 − dˆL )bˆ V (SV (t) + EV (t) + IV (t)) + (1 − dˆL )(1 − m ˆ V )LV (t) ˆ ˆ V LV (t) SV (t + 1) = (1 − dˆV )SV (t)(1 − αˆ V )βR IR (t)/NR (t) + (1 − dˆV )m

ˆ EV (t + 1) = (1 − dˆV )SV (t) 1 − (1 − αˆ V )βR IR (t)/NR (t) + (1 − dˆV )(1 − κˆ V )EV (t)

IV (t + 1) = (1 − dˆV )κˆ V EV (t) + (1 − dˆV )IV (t) SR (t + 1) = SR (t)(1 − αˆ R )βR IV (t)/NR (t) ˆ

498


ˆ IR (t + 1) = (1 − δˆR )SR (t) 1 − (1 − αˆ R )βR IV (t)/NR (t) + (1 − δˆR )(1 − γˆR )IR (t) RR (t + 1) = γˆR IR (t)

ˆ XR (t + 1) = XR (t) + δˆR SR (1 − (1 − αˆ R )βR IV (t)/NR (t) ) + (1 − γˆR )IR (t) . By way of example, we derive in detail the equation for SV (t + 1). First, the expected number of bites made by a susceptible mosquito in a unit time interval is βˆR , and the expected number of times the mosquito bites an infected bird is βˆR IR (t)/NR (t). The probability of a mosquito avoiding the infection arising from a single bite on an infected bird is 1 − αˆ V and hence the probability of a susceptible mosquito avoiding infection in a given time step is (1 − αˆ V )βˆR IR (t)/NR (t) . The number of susceptible mosquitoes remaining susceptible in a given time step, SV (t)(1 − αˆ V )βˆR IR (t)/NR (t) , is augmented by the number of larvae maturing m ˆ V LV (t), and then is diminished by natural mortality which is avoided with probability 1 − dˆV . Using a similar approach, it is possible to derive the other equations above. Here, the parameter constraint for existence of a disease-free equilibrium is bˆ V (1 − dˆV ) = dˆV 1 +

dˆL m ˆ V (1 − dˆL )

.

The disease-free equilibrium is (LV , SV , EV , IV , SR , IR , RR , XR ) =

(1 − dˆL )bˆ V ∗ ∗ ∗ AV , AV , 0, 0, NR , 0, 0, 0 . (1 − dˆL )m ˆ V + dˆL

The infected variables are EV , IV and IR . Linearizing the equations for these variables about the DFE and writing the resulting matrix as F + T, where F includes only new infections and the column sums of T are less than one, gives the following nonnegative matrices: 

A∗ ˆV )βˆR ln(1 − αˆ V ) V 0 0 −(1 − d  NR∗   F = 0 0 0   0 0 −(1 − δˆR )βˆR ln(1 − αˆ R )    T=  

(1 − dˆV )(1 − κˆ V )

0

0

(1 − dˆV )κˆ V

1 − dˆV

0

0

0

(1 − δˆR )(1 − γˆR )

   .  

       


499

Then the matrix F(I − T)−1 is given as

    

0

0

−(1−dˆV )βˆR ln(1−αˆ V )A∗V (δˆR +γˆR (1−δˆR ))NR∗

0

0

0

−(1−δˆR )βˆR ln(1−αˆ R )(1−dˆV )κˆ V −(1−δˆR )βˆR ln(1−αˆ R ) dˆV (dˆV +κˆ V (1−dˆV )) dˆV

    

0

and the spectral radius of F(I − T)−1 is R0 =

(1 − dˆV )2 (1 − δˆR )βˆR2 κˆ V ln(1 − αˆ V ) ln(1 − αˆ R ) A∗V . NR∗ dˆV (dˆV + κˆ V (1 − dˆV ))(δˆR + γˆR (1 − δˆR ))

3.3. Comparison of discrete with continuous time version The discrete-time version in Section 3.2 has a natural correspondence with the continuous time version considered earlier in Section 3.1. In particular, due to the assumed order of events for the discrete system in Section 3.2, the events in step (1) are conditioned upon surviving mortality in the previous step (2). Hence the parameters are conditioned upon surviving natural (dˆV or dˆL ) or disease-induced (δˆR ) mortality. This means that δR is replaced by δˆR , dV is replaced by dˆV , γR is replaced by γˆR (1 − δˆR ), bV is replaced by bˆ V (1 − dˆL ), κV is replaced by κˆ V (1 − dˆV ), and mV is replaced by m ˆ V (1 − dˆL ). The infection rate αV βR is replaced by −(1 − dˆV )βˆR ln(1 − αˆ V ), and infection rate αR βR is replaced by −(1 − δˆR )βˆR ln(1 − αˆ R ). With these substitutions the above basic reproduction number is identical to the previous continuous time version (Section 3.1), and R0 can again be interpreted as a geometric mean. 4. West Nile virus model of Thomas and Urena (2001) 4.1. Discrete-time model We now give the discrete-time model, as formulated by Thomas and Urena (2001), with the time step of 1 week. In the model, the authors include compartments for susceptible, exposed and infectious mosquitoes; susceptible, infectious and recovered birds; and susceptible, infectious and recovered humans. They assume that a proportion of mosquito births is infected, and so goes into the exposed class (vertical WN virus transmission). Birds are assumed to die naturally (not from the virus), human death is omitted, and mass action incidence is assumed. The biting rate parameter, which is defined as ‘the probability that one mosquito bites one bird’ in a given week is constrained to lie between zero and one, unlike the biting parameter of Section 3.2, βˆR , which is defined as the number of bites made per susceptible mosquito per time step. Furthermore, there are no separate transmission

500


parameters αˆ V and αˆ R in the original model (Thomas and Urena, 2001). The probability of virus transmission between mosquitoes and birds is implicitly assumed to be αˆ V = αˆ R = 1. To distinguish the biting rate parameters in this model from the model discussed Section 3.2, we use β˜R and β˜H for the (Thomas and Urena, 2001) model in this section. The order of events is not explicitly stated in (Thomas and Urena, 2001). Control by insecticide spraying every other week is through the introduction of a time dependent function c(t) (death due to spraying), namely c(t) =

0, t even constant c ∈ (0, 1), t odd.

We keep this spraying function in our model formulation, so as to facilitate comparison with the model in Thomas and Urena (2001), although the control level is taken to be equal to zero in the analysis that follows. In the common notation, the dynamics are given by the difference equation system for vectors, reservoirs and humans: Vectors (V): SV (t + 1) = (1 − β˜R ) IR (t) (1 − dˆV + bˆ V − c(t))SV (t) + (1 − ρˆV )bˆ V (EV (t) + IV (t)) EV (t + 1) = (1 − (1 − β˜R ) IR (t) )(1 − dˆV + bˆ V − c(t))SV (t) + (1 − κˆ V )(1 − dˆV − c(t))EV (t) + ρˆV bˆ V (EV (t) + IV (t)) IV (t + 1) = (1 − dˆV − c(t))(IV (t) + κˆ V EV (t)). Reservoirs (R): SR (t + 1) = bˆ R (IR (t) + RR (t)) + (1 − β˜R ) IV (t) (1 − dˆR + bˆ R )SR (t) IR (t + 1) = (1 − (1 − β˜R ) IV (t) )(1 − dˆR + bˆ R )SR (t) + (1 − γˆR )(1 − dˆR )IR (t) RR (t + 1) = (1 − dˆR )RR (t) + γˆR (1 − dˆR )IR (t). Humans (H): SH (t + 1) = bˆ H NH (t) + (1 − β˜H ) IV (t) SH (t) IH (t + 1) = (1 − (1 − β˜H ) IV (t) )SH (t) + (1 − γˆH )IH (t) RH (t + 1) = RH (t) + γˆH IH (t). Note that we have corrected a bracket in their IV (t + 1) equation, and have also omitted a term (1 − κˆ V )bˆ V EV (t) on the right side of the equation for EV (t + 1) as given in Thomas and Urena (2001)(Eq.(3)) as this does not seem biologically reasonable, since vertically transmitted disease is accounted for by the term containing ρˆV . The total number of humans is given by NH = SH + IH + RH .


501

The nonlinearities in the infection terms differ from those assumed in Section 3.2. For example, the nonlinearity in the fraction of susceptible mosquitoes avoiding infection, (1 − β˜R ) IR (t) , differs from the nonlinearity for the discrete model of Section 3.2, (1 − αˆ V )βˆR IR (t)/NR (t) , unless 1 − β˜R = (1 − αˆ V )βˆR /NR (t) , or equivalently βˆR ln(1 − αˆ V )/NR = ln(1 − β˜R ). A similar argument applied to the bird population shows that the nonlinear incidence terms are equal only when βˆR ln(1 − αˆ R )/NR = ln(1 − β˜R ). The model assumptions that lead to the differences in nonlinear transmission terms can be summarized as follows. If a single infected mosquito were introduced into a population of birds, the probability that any given susceptible bird avoids infection during the weekly time-step is 1 − β˜R , according to Thomas and Urena (2001), and is (1 − αˆ V )βˆR /NR (t) , according to Wonham et al. (2004). The first assumes that the probability of avoiding infection is independent of the bird population size, while the second assumes that the probability of avoiding infection is an increasing function of bird population size. Constraints for the existence of a DFE are dˆV = bˆ V , dˆR = bˆ R , bˆ H = 0 and c(t) = 0, and the DFE has SV = NV∗ , SR = NR∗ , SH = NH∗ and all other state variables zero. The model of Thomas and Urena (2001) includes vertical transmission of disease through infected vector births (the parameter ρˆV ). This must also be included in F giving rise to a term in the matrix F. With c(t) = 0 and infected variables EV , IV , IR , IH , this model gives F(I − T)−1 as a 4-by-4 matrix, where the IH variable plays no role in the calculation because humans are dead-end hosts and so do not impact infection of other species (see Appendix for matrices F and T). This matrix is reducible, with two of its eigenvalues being zero, thus its spectral radius is given as the largest eigenvalue of the reduced matrix 

ρˆV

     − ln(1 − β˜R )κˆ V (1 − dˆV )N∗ R dˆV (dˆV + κˆ V (1 − dˆV ))

 − ln(1 − β˜R )NV∗ dR + γˆR (1 − dˆR )   .   0

For ρˆV > 0, some births are infected, then R0 is equal to   ˆV )κˆ V ln2 (1 − β˜R )NV∗ NR∗ (1 − d 1 . ρˆV + ρˆV2 + 4 2 dˆV (dˆV + κˆ V (1 − dˆV ))(dˆR + γˆR (1 − dˆR )) If there are no infected births, then R0 is reduced and is given explicitly as a square root, namely: R0 =

(1 − dˆV )κˆ V ln2 (1 − β˜R )NV∗ NR∗ . dˆV (dˆV + κˆ V (1 − dˆV ))(dˆR + γˆR (1 − dˆR ))

Although the R0 for ρˆV = 0 is superficially different from the R0 for the discretetime version of the Wonham et al. (2004) model in Section 3.2, the two can be

502


closely connected if it is assumed that the value of β˜R can be related to the parameters βˆR , αˆ V , αˆ R , and NR as discussed earlier. Namely, βˆR ln(1 − αˆ V )/NR∗ = ln(1 − β˜R ) and βˆR ln(1 − αˆ R )/NR∗ = ln(1 − β˜R ). Once these assumptions are made, the only difference between the R0 terms is in the precise way that the mortality terms appear in the formula. As demonstrated for the Wonham et al. (2004) model in the Appendix, the manner in which the mortality terms appear in R0 depends upon the precise ordering of events within one-time step of the discretetime model. We believe that this is the reason for the discrepancy between the two R0 terms (for the Thomas and Urena, 2001 model with ρˆV = 0 of this section and for the discrete (Wonham et al., 2004) model of Section 3.2), but have not pursued analysis of this further. This is because, when formulating their model, Thomas and Urena did not precisely specify the ordering of events within a time step.

4.2. Continuous time version of Thomas and Urena (2001) model Assuming a small time step, we replace t + 1 in the original (Thomas and Urena, 2001) model (modified as noted in Section 4.1) with t + t, where t → 0. We replace probabilities by corresponding rates, so that aˆ = at, where a is bV , dV , κV , bR , dR , γR , respectively, and set β˜R = βR t. Then we expand the functions with respect to t using a Taylor series and neglect all higher order terms (such as (t)2 , (t)3 ). The resulting differential equation system, obtained as the limit with t → 0, reads: dSV dt dEV dt dIV dt dSR dt dIR dt dRR dt dSH dt dIH dt dRH dt

= −βR SV IR − dV SV + bV SV + bV (1 − ρV )(EV + IV ) = βR SV IR − dV EV − κV EV + ρV bV (EV + IV ) = −dV IV + κV EV = bR NR − dR SR − βR IV SR = βR IV SR − (dR + γR )IR = γR IR − dR RR = −βH IV SH + bH NH = βH IV SH − γH IH = γH IH .


503

The disease-free equilibrium, with parameter constraints bV = dV , bR = dR and bH = 0, is (SV , EV , IV , SR , IR , RR , SH , IH , RH ) = (NV∗ , 0, 0, NR∗ , 0, 0, NH∗ , 0, 0) . The infected variables are EV , IV , IR and IH but IH can be excluded from the R0 calculations by an argument similar to that used in the original (Thomas and Urena, 2001) model in Section 4.1. Then in variables EV , IV , IR the matrices at the DFE are   ρV bV ρV bV βR NV∗   0 0  F =   0 0 βR NR∗ 0 

dV + κV 0

0

dV

0

 V =  −κV 0

  .

0 dR + γR

Consequently, 

F V −1

    =    

ρV

ρV

0

0

βR κV NR∗ βR NR∗ dV (dV + κV ) dV

 βR NV∗ (dR + γR )      0     0

so the spectral radius of F V −1 is, for ρV = 0, R0 =

βR2 κV NV∗ NR∗ . dV (dV + κV )(dR + γR )

For ρV > 0,   2 ∗ ∗ κ N N β 1 R V V R . R0 = ρV + ρV2 + 4 2 dV (dV + κV )(dR + γR ) With the ratio of A∗V and NR∗ replaced by the product NV∗ NR∗ and δR replaced by dR , the expression for R0 in the continuous (Wonham et al., 2004) model in Section 3.1 is analogous to the above expression with ρV = 0, once we recall that, for the Thomas and Urena (2001) model, αV = αR = 1. As discussed earlier in

504


Section 5, the main difference is due to the different assumptions made about the disease transmission terms. An identification between the R0 expressions for this continuous model and the previous discrete model (Section 4.1) can be made as for the Wonham et al. (2004) models of Section 3. If βR is replaced by − ln(1 − β˜R ), and κV and γR are conditioned on surviving natural death, then the R0 of Section 4.1 is obtained.

5. Discussion and concluding remarks As demonstrated in Section 3, analysis of the change in R0 with model parameters for the Wonham et al. (2004) model predicts that a reduction in mosquito density can be used to control WN virus outbreaks. This prediction is shared by the model of Thomas and Urena for WN virus (2001). However, Wonham et al. (2004) also predicts that a reduction in bird density will actually exacerbate, rather than control, a WN virus outbreak. By way of contrast, this second result is not predicted by the model of Thomas and Urena (2001). In fact, for the Thomas and Urena (2001) model, a reduction in bird density will control, rather than exacerbate, a WN virus outbreak. Clearly, both results cannot be simultaneously correct from a biological perspective. A resolution of this quandary arises from an analysis of how model assumptions shape the formula for R0 . The two models make different assumptions about mosquito-feeding efficiencies, leading to different disease transmission terms. The Wonham et al. (2004) model follows Anderson and May (1991) in using modified mass action terms that assume efficient mosquito searching even when host densities are low (αR βR NSRR IV for mosquitoes to birds and αV βR NIRR SV for birds to mosquitoes). This assumption corresponds to a disease transmission rate that depends only on the proportion of birds susceptible or infected, and is independent of the actual density of birds. Thus a constant disease transmission rate is assumed over a range of bird densities (Fig. 1). The Thomas and Urena (2001) model uses simple mass action terms that assume the encounter rate between mosquitoes and hosts is proportional to host density (in the continuous formulation, these are βR IV SR for mosquitoes to birds and βR SV IR for birds to mosquitoes). This assumption corresponds to a transmission rate that increases linearly with bird density (Fig. 1). As noted earlier, the Thomas and Urena (2001) model does not explicitly state the order of events within in each time step. A further investigation of the Thomas and Urena (2001) model might reformulate the specification of the order of events within each time step, although we do not pursue this here. Since WN-vector mosquitoes in North America typically exhibit a 3-day feeding cycle, it may be reasonable to imagine a saturating functional response of transmission rate to bird density, of which the Wonham et al. (2004) and Thomas and Urena (2001) models each represent a part (Fig. 1). Incorporating such dynamics into a single model would require a different transmission term (McCallum et al., 2001). Although both model assumptions have a sound theoretical basis, they yield starkly different predictions as to the effect of bird control (as calculated from


505

Fig. 1 Two different assumptions about mosquito-feeding efficiency lead to different WN transmission dynamics. The (Wonham et al., 2004) model assumes the transmission rate is invariant across bird densities (horizontal line), whereas the Thomas and Urena (2001) model assumes transmission scales with bird density (diagonal line). When SR = NR = NR∗ , the per mosquito transmission rate to hosts is αR βR in the Wonham et al. (2004) model, whereas the per mosquito transmission rate to hosts is βR NR∗ in the Thomas and Urena (2001) model. Biologically, transmission rates likely exhibit a response that is a combination of the two (solid lines), depending on the equilibrium density of hosts NR∗ .

R0 ) on WN virus. When bird densities are low, the Wonham et al. (2004) model predicts that the remaining birds receive more bites and become local hot spots for disease transmission with each bird having a high probability of becoming infected and passing on the virus. The Thomas and Urena (2001) model, in contrast, predicts that the disease will die out in regions of low bird density. Thus, the R0 of Thomas and Urena (2001) predicts that bird control would be effective in controlling WN, whereas the R0 of Wonham et al. (2004) predicts that it would be counterproductive. Because the R0 involves linearization about the equilibrium SR = NR∗ , the model yielding correct R0 is the one whose functional response of transmission rate to bird density is valid for typical bird densities NR∗ . In other words, if reduction in bird density from NR∗ means no reduction in the overall biting rate, but simply that the remaining birds are bitten more frequently (i.e., NR∗ is in the ‘flat’ region of Fig. 1), then the R0 of Wonham et al. (2004) pertains. If, by way of contrast, reduction in bird density from NR∗ means a concomitant reduction in the overall biting rate (i.e., NR∗ is in the ‘linear’ region of Fig. 1) then the Thomas and Urena (2001) model pertains. This simple example illustrates our primary purpose. That is, to demonstrate clearly how slightly different, but seemingly reasonable assumptions, going into a model formulation for WN virus, can yield very different biological conclusions on the basis of analysis of R0 . Indeed, as we have shown in Section 3 and the Appendix, a change as simple as moving from continuous to discrete-time formulation can yield several plausible discrete-time models, each with a qualitatively different R0 . As such, this paper is designed as a cautionary note which underscores the importance of model formulation and between-model comparison prior to inferring the efficacy of disease management methods. Section 2 and the subsequent explicit calculations address our secondary purpose, namely

506


to set out clearly the calculations for the basic reproduction number for a continuous or discrete-time model, with particular emphasis on models for WN virus. The two models that we have taken from the literature, specifically the continuous time model of Wonham et al. (2004) and the discrete-time model of Thomas and Urena (2001), plus the additional continuous time models by Lord and Day (2001) and by Bowman et al. (2005) are, to our knowledge, the only mathematical models for West Nile virus currently in the literature. We are considering these, a model of St Louis Encephalitis virus by Lord and Day (2001) and other studies on Encephalitis (Kay et al., 1987; Tapaswi and Ghosh, 1999; Tapaswi et al., 1995; Lord and Day, 2001) for further comparisons.

Appendix A.1 Discrete-time version of Wonham et al. (2004) model—version 2 (different assumptions) To see how the ordering of events can affect the model structure, and eventually R0 , we consider a different set of assumptions from that in Section 3.2. The new assumptions are as follows: 1. Disease-induced and natural mortality and birth occur at the beginning of the time step. 2. Infection and transfer occur at the end of the time step. The corresponding model in discrete time is formulated as LV (t + 1) = (1 − m ˆ V )[bˆ V (SV (t) + EV (t) + IV (t)) + (1 − dˆL )LV (t)] βˆR (1−δˆR )IR (t)

SV (t + 1) = (1 − dˆV )SV (t)(1 − αˆ V ) NR (t)−δˆR IR (t) +m ˆ V [(1 − dˆL )LV (t) + bˆ V (SV (t) + EV (t) + IV (t))] βˆR (1−δˆR )IR (t) EV (t + 1) = (1 − dˆV )SV (t) 1 − (1 − αˆ V ) NR (t)−δˆR IR (t) + (1 − dˆV )(1 − κˆ V )EV (t) IV (t + 1) = κˆ V (1 − dˆV )EV (t) + (1 − dˆV )IV (t) βˆR (1−δˆR )IV (t)

SR (t + 1) = SR (t)(1 − αˆ R ) NR (t)−δˆR IR (t) βˆR (1−δˆR )IV (t) IR (t + 1) = (1 − γˆR ) SR (t) 1 − (1 − αˆ R ) NR (t)−δˆR IR (t) + (1 − δˆR )IR (t)

RR (t + 1) = γˆR SR (t) 1 − (1 − αˆ R ) XR (t + 1) = XR (t) + δˆR IR (t).

βˆR (1−δˆR )IV (t) NR (t)−δˆR IR (t)

ˆ + (1 − δR )IR (t)


507

Now, the parameter constraint for existence of a disease-free equilibrium is dˆL bˆ V = dˆV 1 + . m ˆ V (1 − dˆL ) The disease-free equilibrium is (LV , SV , EV , IV , SR , IR , RR , XR ) (1 − dˆL )(1 − m ˆ V )bˆ V ∗ ∗ AV , AV , 0, 0, NR∗ , 0, 0, 0 . = (1 − dˆL )m ˆ V + dˆL The infected variables are EV , IV and IR and  −(1 − dˆV )(1 − δˆR )βˆR ln(1 − αˆ V )

0 0   F = 0 0   0 −(1 − dˆV )(1 − γˆR )βˆR ln(1 − αˆ R )   (1 − dˆV )(1 − κˆ V ) 0 0     . ˆV )κˆ V ˆV T= (1 − d 1 − d 0     0 0 (1 − δˆR )(1 − γˆR )

0

 A∗V NR∗       

0

The matrix F(I − T)−1 is given as  0 0       0 0      −( pˆ V )2 (1 − γˆR )βˆR ln(1 − αˆ R )κˆ V −( pˆ V )(1 − γˆR )βˆR ln(1 − αˆ R ) dˆV (dˆV + κˆ V (1 − dˆV )) dˆV

 pˆ V (1 − δˆR )βˆR ln(1 − αˆ V )A∗V  (δˆR + γˆR (1 − δˆR ))NR∗      0      0

where here pˆ V is written for (1 − dˆV ). Thus the spectral radius of F(I − T)−1 is R0 =

(1 − dˆV )3 (1 − γˆR )(1 − δˆR )βˆR2 κˆ V ln(1 − αˆ V ) ln(1 − αˆ R ) A∗V . NR∗ dˆV (dˆV + κˆ V (1 − dˆV ))(δˆR + γˆR (1 − δˆR ))

Note that this R0 is different than the one calculated for the previous discretetime model in Section 3.2. However, as with the previous model, this R0 can be made identical to the earlier Wonham et al. (2004) R0 , once the parameters are conditioned upon surviving natural (dˆV or dˆL ) or disease-induced (δˆR )

508


mortality, and the infection rates are modified appropriately. For these assumptions, this means that δR is replaced by δˆR , dV is replaced by dˆV , γR is replaced by γˆR (1 − δˆR ), bV is replaced by bˆ V (1 − dˆL ), κV is replaced by κˆ V (1 − dˆV ), and mV is replaced by m ˆ V (1 − dˆL ) as in Section 3.2. The infection rate αV βR is replaced by −(1 − dˆV )(1 − δˆR )βˆR ln(1 − αˆ V ), and infection rate αR βR is replaced by −(1 − δˆR )(1 − γˆR )βˆR ln(1 − αˆ R ). Note that, due to the different set of assumptions regarding the ordering of events, these infection rates are expressed by different parameters than in the previous discrete-time model (Section 3.2): a new factor of 1 − δˆR appears in the αV βR term, and a new factor of 1 − γˆR appears in the αR βR term. A.2 Matrices F and T for Thomas and Urena (2001) model The 4-by-4 matrices F and T used to compute R0 in Section 4.1 are given explicitly as 

ρˆV bˆ V − ln(1 − β˜R )NV∗ ρˆV bˆ V  0 0 0  F = ∗ ˜ 0 − ln(1 − β )N 0 R  R 0 − ln(1 − β˜H )NH∗ 0 

 0 0  0  0

(1 − κˆ V )(1 − dˆV ) 0 0  κˆ (1 − dˆ ) 1 − dˆV 0 V V  T=  0 0 (1 − γˆR )(1 − dˆR ) 0

0

0

0 0 0

   . 

(1 − γˆH )

Acknowledgements We thank two anonymous reviewers for helpful comments that improved the paper. The first author was supported by NSERC Discovery grant, Canada Research Chair and NSERC CRO grant. The second author was supported by PIMS post doctoral fellowship and by KBN grant 2-P03A-002-23. The third author was supported by NSERC Discovery grant and MITACS. The fourth author was supported by NSERC and Honorary Killam postdoctoral fellowships.

References Anderson, R.M., May, R.M. 1991. Infectious Diseases of Humans, Oxford University, Oxford. Bowman, C., Gumel, A.B., van den Driessche, P., Wu, J., Zhu, H., 2005. A mathematical model for assessing control strategies against West Nile virus. Bull. Math. Biol. 67, 1107–1133. Caswell, H., 2001. Matrix Population Models, Sinauer, Sunderland, MA. Diekmann, O., Heesterbeek, J.A.P., 1999. Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation, Wiley, New York.


509

Higgs, S., Schneider, B.S., Vanlandingham, D.L., Klingler, K.A., Gould, E.A., 2005. Nonviremic transmission of West Nile virus, PNAS, 102, 8871–8874. Kay, B.H., Saul, A.J., McCullagh, A., 1987. A mathematical model for the rural amplification of Murray Valley Encephalitis virus in Southern Australia. Am. J. Epidemiol. 125(4), 690–705. Li, C.-K., Schneider, H., 2002. Application of Perron-Frobenius theory to population dynamics. J. Math. Biol. 44, 450–462. Lord, C.C., Day, J.F., 2001. Simulation studies of St. Louis Encephalitis virus in South Florida. Vector Borne Zooton. Dis. 1(4), 299–315. Lord, C.C., Day, J.F., 2001. Simulation studies of St. Louis Encephalitis virus and West Nile viruses: The impact of bird mortality. Vector Borne Zooton. Dis. 1(4), 317–329. McCallum, H., Barlow, N., Hone, J., 2001. How should pathogen transmission be modelled? Trans. Ecol. Evol. 16(6), 295–300. Tapaswi, P.K., Ghosh, A.K., 1999. Dynamics of Japanese Encephalitis—A study in mathematical epidemiology. IMA J. Math. Appl. Med. Biol. 16(10), 1–27. Tapaswi, P.K., Ghosh, A.K., Mukhopadhyay, B.B., 1995. Transmission of Japanese Encephalitis in a 3-population model. Ecol. Modell. 83, 295–309. Thomas, D.M., Urena, B., 2001. A model describing the Evolution of West Nile-like Encephalitis in New York City. Math. Comput. Modell. 34, 771–781. van den Driessche, P., Watmough, J., 2002. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math. Biosci. 180, 29–48. Wonham, M.J., de Camino-Beck, T., Lewis, M.A., 2004. An epidemiological model for West-Nile virus: Invasion analysis and control application. Proc. R. Soc. Lond. B 501–507.