The following is the established format for referencing this article:
Thomas, D., M. Weedermann, L. Billings, J. Hoffacker, and R. A. Washington-Allen. 2009. When to spray: a time-scale calculus approach to controlling the impact of West Nile virus. Ecology and Society14(2): 21. [online] URL: http://www.ecologyandsociety.org/vol14/iss2/art21/Research, part of Special Feature on Catastrophic Thresholds, Perspectives, Definitions, and Applications

When to Spray: a Time-Scale Calculus Approach to Controlling the Impact of West Nile Virus

West Nile Virus (WNV) made its initial appearance in the New York City (NYC) metropolitan area in 1999 and was implicated in cases of human encephalitis and the extensive mortality in crows (Corvus sp.) and other avian species. Mosquitoes were found to be the primary vectors and NYCís current policy on control strategies involved an eradication program that depends on the synchronicity of the summer mosquito populationís increases with the occurrence of cases in humans. The purpose of this paper is to investigate whether this is the most effective control strategy because past mathematical models assumed discrete behavior that is modeled by difference equations for a single summer season was most important to the virusís life cycle. However, both surviving mosquito eggs and surviving migratory birds incubate the virus during the winter, leading to a continuation of infections in the following warmer spring and summer when the birds return and the eggs hatch. Additionally, the virulence of WNV has been observed to fluctuate with changes in temperature toward warmer conditions. Models are required that account for these multi-seasonal dynamics and time-scale calculus is a newly developed method for resolving the behavior of systems that exhibit both discrete and continuous behavior. We found that, although the static states of the new temperature delay model are no different from older models, simulations indicate that the rate of the infection is affected by avian recovery at a lower temperature threshold. Consequently, eradication strategies should consider that controlling mosquitoes during the fall when colder temperatures occur would cause a fast and efficient drop to a disease-free state. This could prove rather more effective than mosquito control in the warmer months.

The West Nile virus (WNV), a member of the genus Flavirus, first made its appearance in the Western Hemisphere, particularly the New York City metropolitan area, in late August and early September of 1999 (Lanciotti et al. 1999, Asnis et al. 2001). The outbreak resulted in 62 reported infections, eight cases of human encephalitis with seven fatalities, and extensive mortality in American crows (Corvus brachrhychos), fish crows (C. osifragus), and several exotic avian species including a Chilean flamingo (Phoenicopterus chilensis) in the Bronx Zoo (Lanciotti et al. 1999, Asnis et al. 2001, LaDeau et al. 2008). At present, WNV can now be found over most of North America, including reports in Canada and Mexico. It was first identified in Uganda in 1937, with sporadic outbreaks in Africa, the Middle East, and eastern Europe through the 1990s (Smithburn 1940, Hubálek and Halouzka 1999). Besides birds, more than 20 000 cases of WNV have been reported in horses in North America (Centers for Disease Control and Prevention (CDC) 2008). Over 62 species of mosquitoes have tested positive for the WNV, but the Genus Culex has been identified as the dominant vector to avian and human hosts (Wonham et al. 2004). LaDeau et al. (2008) have recently elaborated the transmission ecology of WNV (Fig. 1). Mosquitoes usually contract the disease from biting avian carriers in which the virus concentration in the blood has increased substantially (viremia, such as the American crow). Likewise, the reservoir population becomes infected from mosquito bites. Humans and horses can contract the disease, but are considered a dead-end host because they do not produce sufficient virus concentration in the blood to re-infect mosquitoes if bitten (LaDeau et al. 2008). The virus can be transmitted by ingestion of mosquitoes or infected vertebrates (LaDeau et al. 2008).

LaDeau et al. (2008) have noted that researchers need to evaluate further the causes of temporal and spatial heterogeneity in vector abundance, community composition, and their effects on disease dynamics. Additionally, early predictors of thresholds of annual epizootics and human epidemics, including mosquito abundance and winter pathogen survival, are needed, which is the main focus of this paper.

Ever since the emergence of WNV in North America, there have been attempts to mathematically model the dynamics of the virus in order to provide decision makers with control and/or eradication strategies (e.g., Ghosh and Tapaswi 1999, Lord and Day 2001a, 2001b, Thomas and Urena 2001, Choi et al. 2002, Wonham et al. 2004, Bowman et al. 2005, Cruz-Pacheco et al. 2005, Lewis et al. 2006) (see Table 1).

Thomas and Urena (2001) applied a discrete time system to model the interactions between the virus life cycle and the consequent effects on humans. A discrete time system functions by separating time into generations. They determined that the important function to track in the control of WNV is the total number of mosquitoes as a function of time t in weeks (TM(t); Thomas and Urena 2001). Their analysis indicated that the percentage of total mosquitoes eradicated (spray kill rate) had to be larger than a threshold value to result in average mosquito population decline. Spray kill rates below this threshold value resulted only in a temporary decline of mosquito populations, but not a long-term or time-averaged decline (Thomas and Urena 2001). This result allowed policy makers to determine how much spraying would be required to control the virus by estimating the spray kill rate and ensuring it was above the threshold value found in Thomas and Urena (2001).

The Problem

Despite these encouraging results, all previously known models only accounted for dynamics during a single season (Wonham et al. 2006). It has been observed that, during the North American warm season, the rate of mosquito larval development increases, the mosquito population continuously grows in size, and the number of individual feeds of a mosquito increases with increasing temperature and vice versa. Mosquitoes also lay eggs during the warm season, but the hatched populations abruptly end with the first cold snap. The surviving eggs remain dormant over the winter months, leading to a continuation of infections the following warmer spring and summer when the birds return and interactions occur with the hatched eggs. To further complicate this situation, WNV has been observed to become more virulent by replicating itself at a quicker rate during warmer periods and seasons than during cooler periods (Day 2001, Hartman 2002). Consequently, there is a need for subsequent models to capture these multi-seasonal dynamics of WNV that are a combination of both discrete and continuous behavior.

Time-Scale Calculus (TSC) Model

The re-emergence of mosquitoes each year should be modeled over time that is separated into warm and cold seasons. We couple a continuous model of the mosquito population’s dynamics during each summer season with a discrete model that captures winter bird migration to warmer climates and re-emergence of the virus through newly hatched mosquito eggs during the following summer season (Fig. 2). This type of model is called a dynamic equation over a time scale or time-scale calculus (TSC, Hilger 1990).
Difference equations describe population dynamics by considering time-discrete generation periods. A difference equation population model describes the next generation recursively as a function of the previous generation (Thomas and Urena 2001). Difference equations are discrete versions of differential equations, which relate the rate of growth of the population (the derivative) as a function of the population. In the case of differential equations, time is continuous as opposed to discretely labeled generations. Thus, we refer to difference equations as discrete time equations and differential equations as continuous time equations.

With the availability of computing power, difference equation models have increased in use for population modeling. Many of the properties of long-term behavior of solutions for discrete time models have paralleled well-established results from continuous time theory. Instead of proving each parallel result individually, Stefan Hilger introduced the notion of calculus on an arbitrary time scale, which includes continuous time, discrete time, and a mix of both (Hilger 1990). A time scale is simply any closed subset of real numbers with the purpose of developing an equation that evolves over values in this scale. For example, for a differential equation that models population density, the time scale would begin at time equals zero and run over all positive real numbers. In the case of a difference equation model that describes a population of dividing cells, time is discrete. Each time step is the amount of time it takes for a single cell to divide. In this case, the time scale would be positive integers. Time-scale calculus provides a unified theoretical tool for any combination of differential and difference equations. For example, the time scale for a mosquito population would be:

(1)

Bohner and Peterson (2001, 2003) further develop TSC using many of the usual notions of calculus over time scales, including a generalized derivative, a unified set of differentiation rules for finding derivatives (power, product, quotient, and chain rules), and solutions to first-order equations. In fact, general results are continually being developed and extended at the time of the writing of this article.

Consequently, as the infection cycle of West Nile Virus involves a seasonal component, control strategies to manage the disease require consideration of the time scale in Eq. 1. We outlined the use of TSC to understand multi-seasonal effects of WNV to assist policy decisions on the virus in the New York City metropolitan area in a New Scientist article (Spedding 2003). An outcome of this research was a collaborative effort between the policy makers and mathematicians toward the development of a dynamic equation on a time scale that includes both temperate and colder seasons. A mathematical model can evaluate when spraying will be most beneficial for optimal control. For example, the model can identify the times at which the virus’s replication is at a maximum, which is when eradication efforts would have the greatest impact. Therefore, a multi-seasonal model is necessary to evaluate all strategies, not just times during the warm season.

The goal of the dynamic model is to accurately track and predict two of the components of the life cycle of WNV: the mosquito population and the avian host population in a given week. Our model is based on a single core dynamic model developed by Wonham et al. (2006). This core model was the result of a comparison of the construction and resulting epidemiological properties of seven interesting dynamic models (Ghosh and Tapaswi 1999, Lord and Day 2001a, 2001b, Thomas and Urena 2001, Choi et al. 2002, Wonham et al. 2004, Bowman et al. 2005, Cruz-Pacheco et al. 2005, Lewis et al. 2006). Each of these models, including the core model, described the behavior of WNV over a single summer season. We build upon the two dynamic mathematical models based on the core model to capture seasonal dynamics. In the first model, we focus on single season temperature fluctuations and the delay effects of incubating WNV in mosquitoes. The second model attempts to understand the dynamics of West Nile over several seasons.

Time-Scale Model of West Nile

In order to set up a model on a seasonal-dependent time scale, we need to develop a set of time-scale assumptions that identify the beginning and end of each season and the periods when control strategies will be employed. As this is the first attempt to model over a multi-season time scale, we present the simplest case where the mosquitoes and reservoirs are not coupled as the time scale consists of the colder seasons (week 20 to week 52 of the year) when there is no mosquito activity. This creates a jump in time from week 20 (the end of the first warm season) to week 52 (the beginning of the second warm season).

Under these assumptions, the time scale would be: Season 1 = [0, 20], Season 2 = [52, 72], Season 3 = [104, 124] and so forth where time is measured in weeks (Fig. 3). Note that the theoretical mosquito population increases due to dormant eggs that hatch at the beginning of the warm season (Fig. 3). We point out that all the mathematical results on this particular time scale can be easily generalized to any time scale where Season 1 =[0,x], Season 2=[52, x+52], Season 3=[104, x+104] and so forth where x represents the number of weeks in the warmer season. The approach of the unified calculus is to define the rate of change over all seasons by one expression. This expression should reduce to the regular derivative if the seasons were not disconnected and should reduce to discrete change if the time scale consisted of integer moments in time. This generalization is what is known as the delta derivative of a function over a time scale (Bohner and Peterson 2001, 2003). The delta derivative of a function depends on notions from abstract mathematics, so we omit the general definition and refer the interested reader to Bohner and Peterson (2001). However, as it should be, if the time scale is the usual real numbers, that is time runs continuously from zero to infinity, then the delta derivative of a function over the time scale reduces to:

(2)

If the time scale is discrete time 0, 1, 2, ... then the delta derivative of a function over this time scale reduces to the difference of the function at the n + 1st time step and the nth time step or:

(3)

Using the delta derivative we model the total mosquito population density by setting the delta derivative of the total mosquito population [TM(t)] equal to the growth rate of the mosquito population (rM) times the original total mosquito population density[TM(t)] minus the death of the mosquitoes due to pesticide spraying [s times the sum of δ(t,ti)TM(t)]. We make the choice to model the population growth by the simplest possibility first because even direct proportional solutions to time-scale problems have been only recently defined. Likewise, the delta derivative of the total bird population [TB(t)] is equal to the growth rate of the total bird population (rB) times the original bird population:

(4)

The second term on the left-hand side of Eq. 4:

(5)

represents the death of mosquitoes due to pesticide spraying. The parameter, s, is the spray kill rate, which models the percentage of total mosquitoes eradicated due to pesticide spraying. The function, δ(t,ti), is an impulse function that is equal to 1 when spraying occurs at time t = ti, otherwise it is 0. The values of ti can be defined as time values in the time scale when spraying occurs. The total impact to the mosquito population will be the sum of the deaths due to spraying [s times the sum of δ(t,ti)TM(t)]. We remind the reader that the left-hand side of Eq. 4 is the delta derivative of the total mosquito and bird population over the seasonal time scale. The notation is deceptive as the equations appear very simple, yet a closed form solution to this system is still ongoing research.

The first result concerns a no spraying management scenario where the decay or growth of the mosquito population depends on the time scale. The length of the winter months is equal to 32 total weeks. We use the result of the asymptotic behavior of a state variable over a time scale (Gard and Hoffacker 2003) where:

(6)

if this holds, then the total mosquito population will decay to zero. On the other hand, if either the converse set:

(7)

or

(8)

hold, then the mosquito population will grow without bound. This result indicates that control strategies such as pesticide application should be a function of the length of the season.

In order to examine analytically the case where the spray kill rate is not zero, we require the Laplace Transform over a time scale that was defined theoretically by Bohner and Peterson (2001, 2003). The analytical solution for the total mosquito population depends on an exact calculation of the inverse Laplace Transform for our specific equation. Even this simple looking equation on a time scale does not yet have a closed form solution. Although there have been significant advances in Laplace Transform computations, the exact formulation of the inverse Laplace Transform for the non-zero spray kill rate case is still ongoing research (Davis et al. 2007).

Because the seasonal impact on WVN is not a discrete transition but a gradual drop in fall temperature, a closer examination of seasonal temperature fluctuations is required. As stated previously, the WNV incubation length increases as the temperature falls (Day 2001). In order to capture these phenomena, we modeled the virus incubation length using a delay differential equation.

Temperature-Dependent Model

Following the core model in Wonham et al. (2006), we constructed a compartmental susceptible-infected-recovered (SIR) model with a temperature-dependent time lag. The mosquito and bird populations were separated into susceptible, infected, and recovered categories. For example, in order for susceptible birds to become infected, a susceptible bird must be bitten and effectively infected by an infected mosquito. Figure 4 depicts a diagram of the interactions denoted by the variables to be used in the formal mathematical model. In order to derive an appropriate model that includes temperature dependence, we modified the core model by considering a delay in virus incubation periods in the vector as a function of temperature.

Each of the sub-population compartments is called a state variable because we are interested in the state of the sub-population as time changes. For example, the total mosquito population can be divided into the following components: the total mosquito population = the population of mosquitoes susceptible to contracting WNV + the population of mosquitoes exposed to WNV but not yet able to infect a host + the population of mosquitoes that are infective + the mosquitoes who are removed from the system through death. The state variables that appear in the core WNV model represent the susceptible, exposed, and infected mosquitoes (denoted with a subscript V for vector) and the susceptible, infected, and recovered birds (denoted with a subscript R for reservoir) (Table 2).

The core model is a system of interconnected differential equations that correlate the state variables to each other. Each sub-population’s differential equation is developed by the sub-population’s rate of loss and rate of gain. For example, the susceptible mosquito sub-population gains new members through births from each of the total mosquito sub-populations. As a fraction of the births of exposed and infected mosquitoes retain the virus, the gain through birth term is going to be the overall constant birth rate of the total mosquito population times the sum of the susceptible mosquitoes (new births from susceptible mosquitoes) and the proportion of the exposed and infected vectors that do not retain the virus upon birth. Mathematically this is formulated by:

(9)

where bV is the mosquito birth rate and ρV is the fraction of births that retain the infection.

The core model includes an encompassing set of assumptions that include all epidemiological features in past mathematical models for WNV. The assumptions behind the core model are based on observed WNV dynamics involving the vector Culex pipiens (Wonham et al. 2006) and six different reservoir species (Work et al. 1955, Brault et al. 2004, Komar et al. 2005, Langevin et al. 2005, Reisen et al. 1996):

Vector deaths are entirely natural and not disease related.

A fraction of vectors that have contracted the disease transmit the virus to their offspring. The offspring initially is in the exposed class (not effectively infected).

A fraction of infected vectors recover from the virus and have a loss of infectivity.

A fraction of reservoirs have loss of immunity from the virus.

Based on specific reservoir species, the parameters of the model can be chosen to reflect the appropriate scenario. The mean and range of parameter values for six different species were collected from the literature and used for numerical simulations of the core model in (Wonham et al. 2006). For example, in the case of the American crow (C. brachyrhynchos), the probability of surviving infection is 0 and the length of infectivity is 4–6 d (Komar et al. 2005).

In order to include the effects of temperature change in the model, we included a delay term (τ) in the exposed mosquito sub-population. This delay term includes information from a past event that impacts the present situation. Based on the assumption that the incubation of the virus is inversely proportional to temperature (Day 2001), we modify the movement from exposed mosquitoes to infectious mosquitoes by the delay term where τ ranges from 0 to 2 weeks. The resulting model alters the exposed and infectious mosquito equations (see Appendix 1). As all the sub-populations are interconnected, the delay has an impact on all the state variables that can be seen in the stability analysis and numerical simulations.

Mathematical epidemiology uses SIR models to answer several questions. The first is to find the static states, one of which should be disease free. A static state is the state that occurs when there is no change in the state variables as a function of time (the derivative is zero). The static states often are what the dynamic states converge to as time evolves. This fits nicely with the second question: What is the stability and size of the static state? Finally, the speed of convergence to the static states is highly relevant. Reaching a disease-free state in a week as opposed to a decade is an important consideration in designing policy for control of WNV.

In order to answer some of these questions, mathematical epidemiologists compute the basic reproduction number of a disease (Anderson and May 1991). This value is the average number of secondary cases one infected individual will cause. If this number is less than one, the virus will converge to a disease-free state. Conversely, if the basic reproduction number is greater than one, the disease will spread through the population and remain endemic. The size of this value is important because the larger the basic reproduction number, the more severe impact the disease has on the population.

Wonham et al. (2006) found that the basic reproduction number was very sensitive to mosquito and bird population densities. This value was also very sensitive to the type of mosquito–bird interaction. In the modified delay model proposed here, we found that the basic reproduction number was not altered by the time delay caused by temperature drops (Eq. A.2A in the Appendix).

Numerical simulations of the delay model were conducted using MATLAB’s dde23-solver for delay differential equations with constant delays (Shampine and Thompson 2001). Parameters were estimated as the mean values of the known parameter ranges for six different reservoir species, Corvus brachyrhynchos, Turdus migratorius, Cyanocitta cristata, Passer domesticus, Mimus polyglottos, Cardinalis cardinalis, and the vector species Culex pipiens. The parameter values appear in Wonham et al. (2006) and are based on observed field data on WNV with original citations listed in Table 3. Numerical simulations indicated that the rate of convergence to the disease-free equilibrium is a function of the time lag. Figure 5A plots the susceptible, exposed, and infectious classes of mosquitoes with zero time lag and Fig. 5B plots the same quantities for an incubation time lag of 2 weeks. Because the populations are converging to the same disease-free equilibrium, the difference in convergence rates is not dramatically different (Fig. 5).

Table 4 contains data from the simulations for exposed mosquito populations and recovered reservoirs for two different time lags displaying the numerical difference. These numerical results indicate quicker reservoir recovery as the incubation time of the virus in the mosquito lengthens.

The purpose of our model development for the WNV was to include seasonal effects on virus dynamics. If we incorporate longer virus incubation periods as a function of lower temperature, our model shows that the lengthening incubation period does not impact the disease reproduction number, but does affect the rate of convergence to the disease-free equilibrium. Because the virus incubates for a longer time as the temperature drops, we observe faster avian recovery, which speeds up the development of solutions. The results indicate that eradication efforts would be more effective during colder temperature periods relative to warmer periods. In fact, our results suggest that eradication efforts during the fall could prove more effective than in the warmer months.

If we extend the single season model to a multi-seasonal time-scale equation, it is seen that convergence to a disease-free equilibrium is a function of the length of the winter season. A longer winter season extends the range of allowable mosquito growth rate to force a disease-free equilibrium. This implies that less pesticide spraying is required for a longer winter season.

We expect that if a model also includes a slower rate of mosquito larval development and the number of individual feeds of a mosquito decreases when temperature decreases (see Hartman 2004) this convergence will be even quicker. In terms of determining policy decisions regarding West Nile, the impact of this result is to “strike while the iron is hot,” i.e., stronger controls as the temperature drops can have a more profound effect on controlling the development of WNV the following year.

This study successfully applied time-scale calculus to ultimately provide the metropolitan area of New York City with quantitatively derived recommendations on the timing of insect control to eradicate mosquito populations that carried the WNV. It was found that controlling mosquito populations at the end of a season appears to be an effective way to force declines in the virus. Consequently, this study demonstrated the feasibility of using time-scale calculus to study a complex dynamical system where both continuous and discrete behavior of virus, insect, bird, and human populations had to be considered simultaneously. Furthermore, time-scale calculus provides an alternative tool to other approaches for examining dynamic systems such as catastrophe theory and self-organization that are discussed in this Special Feature (e.g., Lockwood and Lockwood 2008).

For this particular work, continued study of the two time-scale models could partially answer many questions about the development of this virus and how to control it. As demonstrated by this study, time-scale calculus models could be used for policy-making decisions, such as testing spraying strategies dependent on a given time scale and time of year that could possibly eradicate the disease. Possible improvements to this model include connecting the delay model with the time-scale model along with bird migration and developing a time-dependent biting rate parameter that slows as the temperature drops. Such a model would require an equation based on two separate time scales, thus implementing a relatively new and open area of research in time-scale calculus.

Responses to this article are invited. If accepted for publication, your response will be hyperlinked to the article. To submit a response, follow this link. To read responses already accepted, follow this link.

This research was funded in part by NSF Biocomplexity Project No. BCS-030846 to Herman H. Shugart, Jr. of the University of Virginia, a Faculty Development Grant from the Office of the Vice-Provost for Faculty Affairs, University of Virginia to R.W-A, and Faculty startup funds from the Department of Ecosystem Science & Management, Texas A&M University and Texas AgriLife Research of the Texas A&M System. Lora Billings was supported by the National Science Foundation under grant DMS-0414087 and the Army Research Office under grant W911NF-06-1-0320.

Asnis, D., R. Conetta, G. Waldman, and A. Teixeira. 2001. The West Nile virus encephalitis outbreak in the United States (1999–2000): from Flushing, New York, to beyond its border. Annals of the New York Academy of Sciences951(1):161–171.

Bugbee, L., and L. Forte. 2004. The discovery of West Nile Virus in overwintering Culex pipiens (Diptera: Culicidae) mosquitoes in Lehigh County, Pennsylvania. Journal of the American Mosquito Control Association20(3):326–327.

Milby, M., and M. Wright. 1976. Survival of house sparrows and house finches in Kern County, California. Bird Banding47(1):119–122.

Oda, T., K. Uchida, A. Mori, M. Mine, Y. Eshita, and K. Kurokawa. 1999. Effects of high temperature on the emergence and survival of adult Culex pipiens molestus and Culex quinquefasciatus in Japan. Journal of the American Mosquito Control Association15(2):153–156.

Sardelis, M., and M. Turell. 2001. Ochlerotatus j. japonicus in Frederick County, Maryland: discovery, distribution, and vector competence for West Nile virus. Journal of the American Mosquito Control Association17(2):137–141.

Walter, N., and C. Hacker. 1974. Variation in life table characteristics among three geographic strains of Culex pipiens quinquefasciatus.Journal of Medical Entomology11(5):541–550.

Wonham, M. J., T. de Camino-Beck, and M. Lewis. 2004. An epidemiological model for West Nile virus: invasion analysis and control applications. Proceedings of the Royal Society and Biological Science271(1538):501–507.