- Email: [email protected]
Modeling dengue outbreaks
Mathematical Biosciences 232 (2011) 87–95
Contents lists available at ScienceDirect
Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs
Modeling dengue outbreaks Marcelo Otero a, Daniel H. Barmak a, Claudio O. Dorso a, Hernán G. Solari a, Mario A. Natiello b,⇑ a b
Departamento de Física, FCEyN-UBA and IFIBA-CONICET, Argentina Centre for Mathematical Sciences, Lund University, Sweden
a r t i c l e
i n f o
Article history: Received 25 September 2010 Received in revised form 22 April 2011 Accepted 25 April 2011 Available online 5 May 2011 Keywords: Epidemiology Dengue Individual based model Compartmental model Stochastic
a b s t r a c t We introduce a dengue model (SEIR) where the human individuals are treated on an individual basis (IBM) while the mosquito population, produced by an independent model, is treated by compartments (SEI). We study the spread of epidemics by the sole action of the mosquito. Exponential, deterministic and experimental distributions for the (human) exposed period are considered in two weather scenarios, one corresponding to temperate climate and the other to tropical climate. Virus circulation, final epidemic size and duration of outbreaks are considered showing that the results present little sensitivity to the statistics followed by the exposed period provided the median of the distributions are in coincidence. Only the time between an introduced (imported) case and the appearance of the first symptomatic secondary case is sensitive to this distribution. We finally show that the IBM model introduced is precisely a realization of a compartmental model, and that at least in this case, the choice between compartmental models or IBM is only a matter of convenience. Ó 2011 Elsevier Inc. All rights reserved.
1. Introduction Dengue fever is a vector-born disease produced by a flavivirus of the family flaviviridae [1]. The main vectors of dengue are Aedes aegypti and Aedes albopictus. The research aimed at producing dengue models for public policy use began with Newton and Reiter [2] who introduced the minimal model for dengue in the form of a set of Ordinary Differential Equations (ODE) for the human population disaggregated in susceptible, exposed, infected and recovered compartments. The mosquito population was not modeled in this early work. A different starting point was taken by Focks et al. [3,4] that began by describing mosquito populations in a computer framework named Dynamic Table Model where later the human population (as well as the disease) was introduced [5]. Newton and Reiter’s model (NR) favours economy of resources and mathematical accessibility, in contrast, Focks’ model emphasize realism. These models represent in dengue two contrasting compromises in the standard trade-off in modeling. A third starting point has been recently added. Otero et al. (OSS) developed a dengue model [6] which includes the evolution of the mosquito population [7,8] and is spatially explicit. This last model is somewhat in between Focks’ and NR as it is formulated as a statedependent Poisson model with exponentially distributed times (a detailed description of our model is given in Section 2).
⇑ Corresponding author. Tel.: +46 46 222 0919; fax: +46 46 222 4010. E-mail address: [email protected] (M.A. Natiello). 0025-5564/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2011.04.006
All modeling approaches have been further developed [9–15]. ODE models have received most of the attention. Some of the works explore: variability of vector population [9], human population [10], the effects of hypothetical vertical transmission of dengue in vectors [11], seasonality [12], age structure [13] as well as incomplete gamma distributions for the incubation and infectious times [16]. Contrasting modeling outcomes with those of real epidemics has shown the need to consider spatial heterogeneity as well [17]. The development of computing technology has made possible to produce Individual Based Models (IBM) for epidemics [18,19]. IBM have been advocated as the most realistic models [19] since their great flexibility allows the modeler to describe disease evolution and human mobility at the individual level. When the results are only to be analysed numerically, IBM are probably the best choice. However, they are frequently presented in a most unfriendly way for mathematicians as they usually lack a formulation (expression in closed formulae) and are – at best – presented as algorithms if not just in words [18]. In contrast, working on the ODE side, it has been possible, for example, to achieve an understanding of the influence of distribution of the infectious period in epidemic modeling [20–22]. IBM have been used to study the time interval between primary and secondary cases [23] which is influenced, in the case of dengue, mainly by the extrinsic (mosquito) and intrinsic (human) incubation period. In this work an IBM model for human population in a dengue epidemic is presented. The model is driven by mosquito populations modeled with spatial heterogeneity with the method introduced in [8] (see Section 2). The IBM model is then used to
88
M. Otero et al. / Mathematical Biosciences 232 (2011) 87–95
examine the actual influence of the distribution of the incubation period comparing the most relevant information produced by dengue models: dependence of the probability of dengue circulation with respect to the mosquito population and the total epidemic size. Exponential, delta (fixed times, deterministic) and experimental [24] distributions are contrasted (Section 3). The infectious period and the extrinsic incubation period is modeled using experimental data and measured transmission rates (human to mosquito) [24]. The IBM model produced is critically discussed. We show that it can be mapped exactly into a stochastic compartmental model of a novel form (see Section 4) thus crossing for the first time the valley separating IBM from compartmental models. This result opens new perspectives which we also discuss in Section 4. Finally, Section 5 presents the conclusions of this work.
2. The model It is currently accepted that the dengue virus does not make any effect to the vector. As such, A. aegypti populations are independent of the presence of the virus. In the present model, mosquito populations are produced by the A. aegypti model [8] with spatial resolution of one block using climatic data tuned to Buenos Aires, a temperate city where dengue circulated in the summer season 2008–2009 [25]. The urban unit of the city is the block (approximately a square of 100 m 100 m). Because of the temperate climate the houses are not open as it is often the case in tropical areas. Mosquitoes usually develop in the center of the block which often presents vegetation and communicates the buildings within the block. The model then assumes that mosquitoes belong to the block and not to the houses and they blood-feed with equal probability in any human resident in the block. A. aegypti is assumed to disperse seeking for places to lay eggs. The mosquito population, number of bites per day, dispersal flights and adult mortality information per block is obtained from the mosquito model [8]. The dynamical evolution of the mosquito population is modeled as a Markov Jump Process. It starts with a list of ‘events’ modifying the populations. Each event occurs randomly according to a Kendall–Feller process [26,27] (approximated here as in [28]) and modifies the subpopulations accordingly (e.g. mature eggs turning to larval state, mosquito bite, oviposition, etc.). Mosquitoes can then perform appentential flights and disperse to lay eggs, in contrast to [5] that considers only one of such urban blocks. The model attempts to use only the minimal number of events that are biologically relevant. The events describe basic biological facts and the values of the associated transition rates and their dependence with e.g., temperature are taken from experimental observation. The stochastic framework is motivated by the intrinsic impredictability of natural processes. An interesting issue in the present model is that it has no adjustable parameters, other than the number of breeding sites (BS) per block. Specifically, the model considers six different populations per block: eggs E, larvae L, pupae P, adult females in their first gonotrophic cycle A1, flyer females F and adult females in subsequent gonotrophic cycles A2. The evolution of the six populations in each block is affected by different possible events: death of eggs, egg development and hatching, death of larvae, pupation, death of pupae, pupal development and adult emergence, completion of gonotrophic cycles, death of young adults A1, oviposition by F flyers, death of A2 adults, death of flyers F and flyers dispersal. We consider that during the gonotrophic cycles the mosquito dispersal is negligible, adults A1 and A2 remain near the breeding sites where they grew, until ready to blood feed. After the blood ingest and the completion of their gonotrophic cycles, adult females begin
to fly (becoming flyers F) in search of oviposition sites. Therefore, every time an adult A1 or an adult A2 becomes a flyer one bite was performed. The rates associated to the Markov jump process connecting the different populations are taken from [7,8]. The model disregards differences in virus strains [1], as well as the clinical complications that arise more often in reincident crossinfections. While multiple virus circulation is possible [29,30] its dynamical relevance has yet to be examinated. In this work we attempt to lay the basis for a flexible model that will later develop further incorporating important dynamical effects such as human’s mobility [31], co-circulation of several viral strains and support the testing of proposed public policies, particularly in view that a vaccine for dengue is not expected in the near future [1]. The model calculates every day the mosquito populations and all the events which occur at rates that depend not only on population values but also on temperature, which in turn is a function of time since it changes over the course of the year. Hence, the dependence on the temperature introduces a time dependence in the event rates. The time-step of the model has been fixed at one day. The human population of each block is fixed in the present work and the disease is spatially spread by the mosquito alone. The evolution of the disease in one individual human, h, proceeds as follows: Day d = d0 The virus is transmitted by the bite of an infected mosquito. The transmission probability is set to pmh = 0.75 [2]. Day d = d0 + sE(h) The human h becomes infective (h is said to be exposed to the virus during this period of time). Day d = d0 + sE(h) + j For 1 6 j 6 sI the human h is infective and transmits the virus to a biting mosquito with a probability phm(j). sI indicates the duration in days of the viremic window, in this article taken to be of 5 days following experimental results in [24]; the values of phm(j) are given in Table 1. Day d > d0 + sE(h) + sI The human h is recovered and no longer transmits dengue. The cycle in the human being is then of the form susceptible, exposed, infected, recovered (SEIR). The virus enters the mosquito when it bites a viremic human with a probability phm(j) depending of the day j in the infectious cycle of the specific human bitten (see Table 1). The cycle continues with the reproduction of the virus within the mosquito (extrinsic period), lasting sm days (in this work sm was set to 8 days). After this reproduction period the mosquito becomes infectious and transmits the virus with a probability pmh when it bites. The mosquito follows a cycle susceptible, exposed, infected (SEI) and does not recover [1,5,24,6]. The adult female mosquito population as produced by the A. aegypti simulation is then split into susceptible, sm stages of exposed and one infective compartment according to their interaction with the viremic human population and the number of days elapsed since acquiring the virus. The epidemic starts when one or more humans become viremic. The algorithm followed is: 1. Give individual attributes, sE(h) according to prescribed distribution. 2. Read geometry, size of viremic window, probabilities phm(j), j = 1, . . . , sI, human population in each block, day of the year when the epidemic starts. Table 1 Human-to-mosquito dengue virus transmission probabilities for the different days of the infected stage.
Probability
Day 1
Day 2
Day 3
Day 4
Day 5
0.25
0.8
1
0.95
0.5
89
M. Otero et al. / Mathematical Biosciences 232 (2011) 87–95
6. Repeat all over again from (5a). Each iteration is a new day of the simulation. We take the opportunity to observe that the differences between [6] and the present model are more in the algorithmic side than in the biological content. The present IBM code is substantially more efficient as well as flexible than the compartmental model in [6], allowing for larger cities and better statistics. On
the biological side, the present model is more realistic as it incorporates the experimental values for phm (exponentially distributed in [6]) as well as a more realistic distribution of the extrinsic incubation time.
0.5
Relative frequency
0.4
0.3
0.2
0.1
0
0
1
2
3 4 5 Intrinsic Incubation Period / days
6
7
8
0.4
P
0.6
0.8
1.0
Fig. 1. Distribution for the duration of the exposed period, according to Nishiura’s data [24].
0.0
0.2
N E1 E2 D
0
200
400
600
800
1000
800
1000
0.4
P
0.6
0.8
1.0
BS
0.2
N E1 E2 D
0.0
3. Initialise the blocks with the human population. 4. Read total adult female mosquito population of the day (M), bites, flights to neighbouring blocks and mosquito death probability. Initialise all mosquitoes as susceptible (MS). Set infective bites to zero. 5. Day-loop begins: (a) Calculate the amount of surviving infected mosquitoes. Compute surviving mosquito population with d exposed days, age them by one day. Exposed mosquitoes evolve to infected ones after d = sm days (Use binomial random number generator). (b) Compute probability for a mosquito bite to transmit dengue as pminf = pmhMI/M (compound probability of being infective and being effective in the transmission). Each bite is an independent event according to the underlying mosquito model. (c) Compute the probability for a bite to be made by a susceptible mosquito out of all the non-infective bites potras = MS/ (M(1 pminf)). The non-infective bites have probability (1 pminf), and come from infected mosquitoes failing to transmit the virus, exposed and susceptible mosquitoes, M (1 pminf) = M pmhMI = MS + ME + (1 pmh)MI. (d) Calculate the number of infected and susceptible mosquito bites using binomials with the previous probabilities and the number of total bites. (e) Compute number of humans bitten in each day of the infected state, HI(j). (f) Compute the number of new exposed mosquitoes, taking into account that the probability of human-mosquito contagion is dependent on the stage-day of the human infection. The amount of new infected insects is chosen using binomials. For this purpose, we add the results of the calculation of the number of susceptible mosquitoes that bite humans and PI get infected, M NE ¼ sj¼1 Binðphm ðjÞ; HI ðjÞÞ. Where MNE are the new exposed mosquitoes, Bin(HI(j), phm(j)) is a binomial realization with the day-dependent probability phm(j) and HI(j) the quantity of infected humans bitten by susceptible mosquitoes on infection day j. (g) Perform a random equi-distributed selection of humans bitten by infectious mosquitoes. Build a table of bitten individuals. (h) Update the state of all the humans. If the human belongs to the Susceptible state and has been bitten according to the table built in (5 g) then change the state of selected human to Exposed, record d0 for each exposed human. Susceptible individuals not bitten by an infective mosquito will remain as such and consequently their intrinsic time d remains in 0. Increase the intrinsic time of Exposed and Infected humans by one. Those Exposed individuals for which their intrinsic time has surpassed the value d = d0 + sE(h) are moved to the Infected state, while those Infected individuals whose intrinsic time is larger than d = d0 + sE(h) + sI, are moved to the Removed state (i) Compute the number of individuals in each state for every cell. (ii) Read the total adult female mosquito population of the day (M), bites, flights to neighboring blocks.
0
200
400
600
BS Fig. 2. Probability of local circulation of virus (probability of having at least one secondary case). Top: with a tropical temperature scenario, Bottom: with a temperate climate, BS: Number of breeding sites per block.
90
M. Otero et al. / Mathematical Biosciences 232 (2011) 87–95
3. Epidemic dependence on the distribution of the exposed period We implemented four different distributions for the duration of the exposed period assigned to human individuals: Nishiura’s experimental distribution [24], a delta and exponential distributions with the same mean that the experimental one and an exponential distribution with the same median than the experimental one. We call them N, D, E1 and E2, respectively. See Fig. 1 for an histogram of Nishiura’s experimental distribution. The study was performed in two different climatic scenarios, one with constant temperature of 23 degrees Celsius, that represents tropical regions and one with the mean and amplitude characteristic of Buenos Aires, a city with temperate climate. The number of effective breeding sites [7] was varied between 50 and 1000. The simulations were performed using identical mosquito populations in all cases (same data file), for an homogeneous urban area of 20 by 20 blocks, hosting 100 people per block, making a total of 40,000 human beings. Hence, all differences correspond to the disease dynamics that was previously identified as the main source of stochastic variations. The simulations with temperate climate were started on January 1st, i.e., ten days after the summer solstice (December 21st in the southern hemisphere). The reported statistics is computed by averaging the outcomes of 1000 simulations with different seeds for the pseudo-random routines. The mosquito part of the simulation has been successfully tested in larger modeling environments, suggesting that just mosquito dispersal is an important ingredient contributing to the mosquito
persistence (survival) in cities with temperate climate [8]. The human population in this model is taken to be static, with fixed population density, hence a larger population size will not change the overall results. The effects of moderately changing the human population can be easily understood. The window for secondary cases will not change (does not depend on the population), the epidemic size will be proportional to the human population size, while the probability of virus-circulation as well as R0 decrease when the population increases (the probability of the index case to be bitten by a mosquito varies as the inverse of the population size [32]), we therefore do not present confirmatory results of these well-understood effects. The main questions were: considering the total number of recovered individuals, how is the distribution of epidemic sizes influenced by the choice of distribution? How does the probability for having no secondary cases change? How is the predicted duration of the epidemic influenced by our choices? And finally, how is the distribution of time between epidemiologically related cases affected? This time is called the generation time in [33]. The epidemiologically relevant question is: can a given (dengue) case be epidemiologically related to (one of) the previously detected case(s)? The proper answer to the later question will be of the type: considering the time interval separating both events, the probability of they being epidemiologically related is . . . Generation times need to be studied at the probability distribution level. We will present the minimal information of the distributions in the form of box-plots, showing the median, lower and upper quartiles, and the ‘outlier points’. Before showing our results, it is worth to realise that the probability of having no secondary cases will not be sensitive
E1
30000
30000
30000
0 00
250
0 200
10000
150
10000
100
20000
50
20000
10
Epidemics Size
30000
50
250
E2 40000
10
250
D
Epidemics Size
200
BS
40000
BS
200
BS
150
10
150
250
0 200
0 100
10000
50
10000
150
20000
100
20000
50
Epidemics Size
40000
10
Epidemics Size
N 40000
BS
Fig. 3. Box plot graphs for the epidemic size (total number of infected humans) at constant temperature. Top-left: N-distribution, top-right: E1-distribution, bottom-left Ddistribution and bottom-right: E2-distribution. BS: Number of breeding sites per block.
91
M. Otero et al. / Mathematical Biosciences 232 (2011) 87–95
to the choice of distribution for the exposed period, as this probability depends only on the probability of the introduced case being bitten by the mosquitoes, the probability of the mosquitoes of acquiring the virus, surviving the extrinsic period and finally transmitting the virus to a human in a bite. Nothing in this process depends on the choice of distribution. In contrast, we expect the distribution of times between epidemiologically related cases to depend strongly on the choice of distribution, since it reflects the sum of the two incubation periods (intrinsic and extrinsic). Some intuition of the epidemiological situation described can be achieved calculating R0, the basic reproductive number. We follow the definition in [32] and define R0 as the average number of secondary cases produced by a single case when the epidemic starts. A rapidly growing epidemy will have R0 > 1 while a disease going to extinction will display R0 < 1. A plot of this R0 as a function of the number of breeding sites is displayed in Fig. 5. In the two conditions studied in this manuscript (summer conditions in temperate climate and a year-long favourable tropical climate), R0 > 1 at values of the breeding sites per block as low as 15. Confirming the reasoning above, there is no sensitivity for the probability of having local circulation of the virus (defined as having at least one secondary case following the introduced case) with the statistical distribution of the exposed period, see Fig. 2. However, we note that our model does not force local circulation of the virus, since the probability of circulation is not identically equal to one. According to the Law of Large Numbers, in a sufficiently large number of simulations, a portion of them will yield extinction of the disease.
The size of the epidemic at constant temperatures makes a transition from very small outbreaks to a large outbreak reaching almost all the people in the simulation. The transition happens in the region 50–100 breading sites per block for the four distributions studied, see Fig. 3. We detected no epidemiologically important differences produced by the use of one or another distribution function. The size of epidemics with seasonal dependence is presented in Fig. 4. The results for the D-distributions and N-distribution do not present differences. The E1-distribution (equal mean) overestimates the final size while the E2-distribution substantially agrees with the N (experimental) one. Both exponential distributions present a larger variance than the N-distribution. This difference may matter when worst-possible scenarios are considered. The size of the epidemic outbreak begins to increase with the number of breeding sites in the 100–150 region for all the distributions. It is possible to argue that the differences in epidemic size observed for different distributions in the case with seasonal dependence correspond to a faster evolution of the epidemic outbreak during the time-window of favourable conditions. On the other hand, in the constant temperature scenario the epidemic outbreak stops because of the decrease in susceptible people produced by the epidemic. In this case we observe no significant differences, see Fig. 6, the epidemic size always ends up around 40,000 individuals. However, the epidemics evolves faster the larger the number of breeding sites (rightmost plots). It is worth observing that major outbreaks last more than one year in this 40,000 people urbanization. Since there is no human movement incorporated, the dura-
25000
BS
900
950
1000
950
1000
850 850
900
750
800 800
700 700
750
600
650
600
650
500
550
500
450
400
350
300
250
200
950
1000
900
850
800
750
700
650
600
550
500
450
0 400
0 350
5000
300
5000
150
10000
100
10000
15000
50
15000
20000
10
Epidemics size
20000
550
E2
25000
250
450
D
25000
200
350
BS
30000
150
400
BS
30000
100
250
10
950
1000
850
900
800
750
650
700
550
600
500
400
450
300
350
200
250
150
0 100
0 10
5000
50
5000
300
10000
200
10000
15000
150
15000
20000
50
20000
100
Epidemics size
25000
10
Epidemics size
E1 30000
50
Epidemics size
N 30000
BS
Fig. 4. Box plot graphs for the epidemic size with seasonal dependence. Top-left, N-distribution, top-right, E1-distribution, bottom-left D-distribution and bottom-right, E2distribution. BS: Number of breeding sites per block.
92
M. Otero et al. / Mathematical Biosciences 232 (2011) 87–95
30
tious material, and the appearance of symptoms. . . Several implications can be drawn from these epidemiological principles and serve as a conceptual model for an individual infection process. First, it is most appropriate to represent the individual infection as a series of discrete events and periods. Second, the discrete events have no duration in their own right and only trigger the change between infection periods. Third, the discrete periods indicate the infection status of an individual and are part of the individual’s characteristics.
D N E1 E2
25
R0
20
15
10
5
0
0
200
400
600
800
1000
800
1000
BS
30
D N E1 E2
25
R0
20
15
10
5
0
0
200
400
600
BS Fig. 5. Estimation of the basic reproductive number R0 as a function of the number of breeding sites per block. Top: with a tropical temperature scenario. Bottom: with a temperate climate.
tion of the outbreak depends critically on the dispersion of female mosquitoes. At this point we could ask: is there any relevant statistics that depends on the distribution of exposed times? The answer is yes. Assume a case of imported dengue is detected, how long do we have to wait to know if we are facing an outbreak or not? The time elapsed between the primary and the first secondary case is given by adding the extrinsic incubation period and the exposed time. Box plots produced after 1000 simulations of outbreaks wit a mosquito population supported on 200 BS per block are displayed in Fig. 7. We observe that in this case the exponential distributions exaggerate the dispersion of results producing too early as well as too late cases as compared with the experimental distribution while the delta-distribution compresses the ‘alert window’ too much with respect to the experimental distribution. 4. Is the IBM model an implementation of a compartmental model? The use of IBM in epidemiology is usually advocated on an ontological perspective [34,19], we quote the argument in [19] . . .the epidemiology literature has always described an infection history as a sequence of distinct periods, each of which begins and ends with a discrete event. The critical periods include the latent, the infectious, and the incubation periods. The critical events include the receipt of infection, the emission of infec-
The evolution of dengue at the individual level is described in detail in the literature, including experimental results [24] which are seldom available for other illnesses. Thus, dengue is a good case to put the thesis at test. The description in terms of events immediately calls for stochastic population models, while the different human to mosquito transmission probabilities can be handled easily introducing age structure in the infective human population, The description of the proper probability distribution for the exposed (latent) period is the major obstacle towards a compartmental stochastic model. Each human individual spends k days, k 2 {1, . . . , sE}, in the exposed state before becoming infective. The probability of any individual to spend k days is P(k). Hence, in the group of N individuals that become infected at d = d0, a portion N(k) will spend k days before they become infective, where N(k) is a random deviate taken from Multinomial(N, P(1), . . . , P(sE)), a multinomial distribution. This is the only information available to the simulation. IBM generate additional structure, since each individual is assigned to one of the classes N(1), . . . , N(sE) in an entirely random way. This apparent additional information is in fact arbitrary because of the randomness and it is averaged out when presenting the results. Although the algorithm of the IBM model assigns the time spent in the exposed class to the individual, the final cause of having a distribution of exposed times is not known (neither to us nor to the algorithm). Different exposed times may arise either because of differences (in virus resistance) in the exposed individuals or because of differences in the incoming virus due e.g., to biological processes concerning the development and transmission of the virus by the mosquito. Let us define
b PðkÞ ¼
sE X
PðjÞ;
ð1Þ
j¼k
QðkÞ ¼
PðkÞ ¼ Pðj ¼ k=j P kÞ; b PðkÞ
bk ¼ N N
k X
NðjÞ:
ð2Þ
ð3Þ
j¼1
Clearly Q(k) 6 1 and Q(sE) = 1. Consider the numbers E(1) = N, E(k) = Binomial(E(k 1), 1 Q(k)) and N(k) = E(k 1) E(k) for 2 6 k 6 sE. We will now recall a few elementary results regarding multinomial distributions. We write Multi for the multinomial distribution. Lemma 1. Multinomial splitting
MultiðN; pð1Þ; . . . ; pðk 1Þ; pðkÞ; . . . ; pðsE ÞÞ ¼ b ¼ Multi N; pð1Þ; . . . ; pðk 1Þ; PðkÞ b Multi NðkÞ; Q ðkÞ; . . . ; Q ðsE Þ :
ð4Þ
Proof. The proof is simple algebra simplifying the expression for the probabilities in the right side of the equation. h
93
M. Otero et al. / Mathematical Biosciences 232 (2011) 87–95
0 500
1500
2500
3500
0
200
400
600
800
0
200
400
600
800
0
200
400
600
800
0 500
1500
2500
3500
0
200
400
600
800
0
200
400
600
800
0
200
400
600
800
0 500
1500
2500
3500
0
200
400
600
800
0
200
400
600
800
0
200
400
600
800
0 500
1500
2500
3500
0
200
400
600
800
0
200
400
600
800
0
200
400
600
800
Fig. 6. Histogram of the duration (in days) of the epidemic outbreak at constant temperature for different number of breeding sites, left to right 50, 100, 150 and 200 breeding sites (BS) per block. Top line, N-statistics, second line E1-statistics, third line E2-statistics and bottom line D-statistics. x-axis: Duration (in days); y-axis: Relative frequency.
Corollary 1 (Binomial decomposition). The multinomial distribution can be fully decomposed in terms of binomial distributions in the form
MultiðN; pð1Þ; . . . ; pðsE ÞÞ ¼
sE Y
BinomialðEðkÞ; QðkÞÞ:
ð5Þ
k¼1
Proof. Repeated application of the previous lemma is all what is needed. h We now state the application of these results to our modeling problem. Theorem 1. Compartmental presentation Consider sE compartments associated to the days k = 1, . . . , sE after receiving the virus from the mosquito for those humans that are not yet infective. Let the population number in the k compartment be E(k), and the number of humans that become infective on day k be N(k) which is a random deviate distributed with Binomial(E(k), Q(k)), with the definitions given above. Then the exposed period that corresponds to the individuals is distributed with P(k).
Proof. It follows immediately from the corollary.
h
The resulting scheme for the compartmental model is illustrated in Fig. 8. 5. Summary, discussion and conclusions We have developed an IBM model for the evolution of dengue outbreaks that takes information from mosquito populations sim-
ulated with an A. aegypti model and builds thereafter the epidemic part of the evolution. The model was used to explore the actual influence of the distribution of exposed time for humans in those characteristics of epidemic outbreaks that matter the most: determining the level of mosquito abundance that makes unlikely the occurrence of a dengue outbreak and determining the size and time-lapse of the outbreak. The distributions used are (a) an experimentally obtained distribution (Nishiura [24]) (labeled N), (b) and (c) exponential distributions adjusted to have the same mean, respectively the same median as N, labeled E1 and E2, and (d) a fixed-time distribution equal to the experimental mean (labeled the D-distribution). 5.1. Main findings For mosquito populations sufficiently large to support the epidemic spread of dengue, most of the stochastic variability is provided by the epidemic process rather than by stochastic fluctuations in the mosquito population. Thus, the splitting of the models improves substantially the performance of the codes. The probability of producing one or more secondary cases after the arrival of an infective human does not depend on the choice of distribution. The characteristic size of the epidemics under a temperate climate are exaggerated by the E1-distribution but presents no substantial difference for the other distributions. The dispersion of values is exaggerated by both exponential distributions. We observed no important differences in the duration of epidemics developed under a constant temperature since the outbreaks reach almost all the population and the velocity is regulated by the dispersion of the mosquitoes in the absence of movement by humans.
N
E1
60 30 20 10 0
0
0
10
10
20
20
30
30
40
40
50
50
60 50 40
50 40 30 0
10
20
time / days
60
M. Otero et al. / Mathematical Biosciences 232 (2011) 87–95
60
94
E2
D
Fig. 7. Distribution of times between the arrival of an infectious person and the time for the first secondary case. The scenario corresponds to 200 breeding sites per block at constant temperature of 23 °C and 100 people per block. Results correspond to the delta-distribution (D), the experimental distribution (N), the exponential E2 and E1 distributions.
S
E(1)
1−Q(1)
Q(1)
E(2)
E(Te)
Q(Te)=1
Q(2)
Pmh I(1)
1
I(2)
I(Ti)
1 R
Phm(1) 1 IM
EM(Tm)
EM(1)
SM
Phm(Ti)
Fig. 8. Scheme for the progression of dengue. Circles for human subpopulations (susceptible, exposed by day, infective by day and recovered) and squares for mosquito subpopulations. Solid directed arrows indicate daily progression, dotted directed arrows indicate several days of progression, bi-directed arrows with a circle indicate interactions resulting in transitions for members of one population. Probabilities are indicated next to arrows. The mosquito dynamics (birth, death, . . . ) is not represented.
The only statistic able to discriminate easily between the four distributions of exposed time was found to be the time of appearance of the first secondary case (generation time). A result that was expected as well. This fact allows to conjecture that in an extreme situation for the mosquito development, e.g., approaching winter temperatures in a temperate climate, the different delays given by the four distribution will discriminate different behaviours. In the limit, if the incubation time for the first infected case were so long that there are no mosquitos around to pick up the virus, then no disease circulation will occur. That could only happen for extremely late (seasonally speaking) primary infections. 5.2. Limitations The split between mosquito evolution and epidemic evolution is not perfect since the events bite and flight are treated as independent events while they are in fact correlated [35–38], both being related to oviposition.
A general feature of stochastic models is the fact that populations are subdivided into compartments. Within each class, all individuals behave as ‘average individuals’. Extensive model-testing and falsation attempts are the only tool to gain confidence in the proposed structure.
5.3. Concluding remarks The model we developed is shown to be an IBM implementation of a compartmental model, of a form not usually considered. At the level of the description in this work, IBM does not play a fundamental role, contrary to what it has been previously argued [34,19]. Rather, its use is a matter of algorithmic convenience. This result indicates that ‘IBM versus compartmental models’ is not a fundamental dichotomy but it may be a matter of choice (depending of the skills and goals of the user): IBM facilitate coding, compartmental models lend themselves to richer forms of analysis.
M. Otero et al. / Mathematical Biosciences 232 (2011) 87–95
In conclusion, only very specific matters seem to depend on the characteristics of the distribution of exposed times for human beings. The question: do the present results support the inclusion of highly detailed incubation periods? cannot be answered without further specification. The statistical detail required in a given simulation depends on the questions that we intend to answer with the model and on the use to be given to these answers. In the present work the answer to the question: what could go wrong if we omit (or misrepresent) the distribution of exposed times? has been explored. For example, an accurate description of the distribution of exposed times is mandatory to assess if a (supposedly secondary) dengue case is epidemiologically related to a previous one. Although other features of the dynamics appear to be less sensitive to this distribution, we cannot expect this always to be the case. Every omission represents a potential liability but on the other hand, in any conceivable model, detail is always finite. Perhaps one of the most striking and unexpected results of the present work is that the exponential distribution with the same median as the experimental data is to be preferred in front of the one with the same mean, contrary to the popular belief (which was indeed what we believed when we started the work). Looking towards the past, conclusions reached using exponential and delta distributions cannot be objected on such basis. Looking towards the future, simple compartmental models can be constructed as well using realistic distributions and there is no reason to limit the models to the choice: exponential, gamma or delta. References [1] D.J. Gubler, Dengue and dengue hemorrhagic fever, Clin. Microbiol. Rev. 11 (1998) 480. [2] E.A.C. Newton, P. Reiter, A model of the transmission of dengue fever with an evaluation of the impact of ultra-low volume (ulv) insecticide applications on dengue epidemics, Am. J. Trop. Med. Hyg. 47 (1992) 709. [3] D.A. Focks, D.C. Haile, E. Daniels, G.A. Moun, Dynamics life table model for Aedes aegypti: Analysis of the literature and model development, J. Med. Entomol. 30 (1993) 1003. [4] D.A. Focks, D.C. Haile, E. Daniels, G.A. Mount, Dynamic life table model for Aedes aegypti: Simulations results, J. Med. Entomol. 30 (1993) 1019. [5] D.A. Focks, D.C. Haile, E. Daniels, D. Keesling, A simulation model of the epidemiology of urban dengue fever: Literature analysis, model development, preliminary validation and samples of simulation results, Am. J. Trop. Med. Hyg. 53 (1995) 489. [6] M. Otero, H.G. Solari, Mathematical model of dengue disease transmission by Aedes aegypti mosquito, Math. Biosci. 223 (2010) 32. [7] M. Otero, H.G. Solari, N. Schweigmann, A stochastic population dynamic model for Aedes aegypti: Formulation and application to a city with temperate climate, Bull. Math. Biol. 68 (2006) 1945. [8] M. Otero, N. Schweigmann, H.G. Solari, A stochastic spatial dynamical model for Aedes aegypti, Bull. Math. Biol. 70 (2008) 1297. [9] L. Esteva, C. Vargas, Analysis of a dengue disease transmission model, Math. Biosci. 150 (1998) 131. [10] L. Esteva, C. Vargas, A model for dengue disease with variable human population, J. Math. Biol. 38 (1999) 220. [11] L. Esteva, C. Vargas, Influence of vertical and mechanical transmission on the dynamics of dengue disease, Math. Biosci. 167 (2000) 51. [12] L.M. Bartley, C.A. Donnelly, G.P. Garnett, The seasonal pattern of dengue in endemic areas: Mathematical models of mechanisms, Trans. R. Soc. Trop. Med. Hyg. 96 (2002) 387. [13] P. Pongsumpun, I.M. Tang, Transmission of dengue hemorrhagic fever in an age structured population, Math. Comput. Model. 37 (2003) 949.
95
[14] K. Magori, M. Legros, M.E. Puente, D.A. Focks, T.W. Scott, A.L. Lloyd, F. Gould, Skeeter buster: A stochastic, spatially explicit modeling tool for studying Aedes aegypti population replacement and population suppression strategies, PLoS Negl. Trop. Dis. 3 (9) (2009) e508, doi:10.1371/journal.pntd.0000508.. [15] M.L. Fernández, M. Otero, H.G. Solari, N. Schweigmann, Eco-epidemiological modelling of Aedes aegypti transmitted diseases. case study: Yellow fever in Buenos Aires 1870–1871, Preprint available from the authors, submitted for publication. [16] G. Chowella, P. Diaz-Dueñas, J. Miller, A. Alcazar-Velazco, J. Hyman, P. Fenimore, C. Castillo-Chavez, Estimation of the reproduction number of dengue fever from spatial epidemic, Math. Biosci. 208 (2007) 571. [17] C. Favier, D. Schmit, C.D.M. Müller-Graf, B. Cazelles, N. Degallier, B. Mondet, M.A. Dubois, Influence of spatial heterogeneity on an emerging infectious disease: The case of dengue epidemics, Proc. R. Soc. (London): Biol. Sci. 272 (1568) (2005) 1171. [18] V. Grimm, Ten years of individual-based modelling in ecology: What have we learned and what could we learn in the future?, Ecol Model. 115 (2–3) (1999) 129. [19] L. Bian, A conceptual framework for an individual-based spatially explicit epidemiological model, Environ. Plann. B: Plann. Des. 31 (3) (2004) 381. [20] A. Lloyd, Destabilization of epidemic models with the inclusion of realistic distributions of infectious periods, Proc. R. Soci. London B 268 (2001) 985. [21] A.L. Lloyd, Realistic distributions of infectious periods in epidemic models: Changing patterns of persistence and dynamics, Theor. Popul. Biol. 60 (1) (2001) 59. [22] Z. Feng, D. Xu, H. Zhao, Epidemiological models with non-exponentially distributed disease stages and applications to disease control, Bull. Math. Biol. 69 (2007) 1511. [23] P. Fine, The interval between successive cases of an infectious disease, Amer. J. Epidemiol. 158 (11) (2003) 1039. [24] H. Nishiura, S.B. Halstead, Natural history of dengue virus (denv)-1 and denv-4 infections: Reanalysis of classic studies, J. Infect. Dis. 195 (2007) 1007. [25] A. Seijo, Y. Romer, M. Espinosa, J. Monroig, S. Giamperetti, D. Ameri, L. Antonelli, Brote de dengue autóctono en el area metropolitana Buenos Aires. experiencia del Hospital de enfermedades infecciosas F.J. Muñiz, Medicina, 0025-7680 69 (2009) 593. [26] D.G. Kendall, An artificial realization of a simple birth-and-death process, J. R. Stat. Soc., Ser. B (Methodological) 12 (1950) 116. [27] W. Feller, On the integro-differential equations of purely discontinuous Markoff processes, Trans. Am. Math. Society 48 (3) (1940) 488. http:// www.jstor.org/stable/1990095. [28] H.G. Solari, M.A. Natiello, Stochastic population dynamics: The Poisson approximation, Phys. Rev. E 67 (2003) 031918. [29] R.M. Myers, M.J. Varkey, R. Reuben, E.S. Jesudass, Dengue outbreak in Vellore, Southern India, in 1968, with isolation of four dengue types from man and mosquitoes, Indian J. Med. Res. 58 (1) (1970) 24. [30] R.M. Myers, M.J. Varkey, R. Reuben, E.S. Jesudass, B. Benjamin, The 1968 outbreak of dengue in Vellore Southern India, Am. J. Public Health 61 (7) (1971) 1379. [31] D.H. Barmak, C.O. Dorso, M. Otero, H.G. Solari, Dengue epidemics and human mobility, arXiv:1102.3869v1[q-bio.PE]. [32] A.L. Lloyd, J. Zhang, A.M. Root, Stochasticity and heterogeneity in host vector models, Interface 4 (2007) 851. [33] A. Svensson, A note on generation times in epidemic models, Math. Biosci. 208 (1) (2007) 300. [34] J.S. Koopman, J.W. Lynch, Individual causal models and population system models in epidemiology, Am. J. Public Health 89 (8) (1999) 1170. [35] M. Wolfinsohn, R. Galun, A method for determining the flight range of Aedes aegypti (linn.), Bull. Res. Counc. Israel 2 (1953) 433. [36] P. Reiter, M.A. Amador, R.A. Anderson, G.G. Clark, Short report: Dispersal of Aedes aegypti in an urban area after blood feeding as demonstrated by rubidium-marked eggs, Am. J. Trop. Med. Hyg. 52 (1995) 177. [37] L.E. Muir, B.H. Kay, Aedes aegypti survival and dispersal estimated by markrelease-recapture in Northern Australia, Am. J. Trop. Med. Hyg. 58 (1998) 277. [38] J.D. Edman, T.W. Scott, A. Costero, A.C. Morrison, L.C. Harrington, G.G. Clark, Aedes aegypti (diptera culicidae) movement influenced by availability of oviposition sites, J. Med. Entomol. 35 (4) (1998) 578.
Contents lists available at ScienceDirect
Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs
Modeling dengue outbreaks Marcelo Otero a, Daniel H. Barmak a, Claudio O. Dorso a, Hernán G. Solari a, Mario A. Natiello b,⇑ a b
Departamento de Física, FCEyN-UBA and IFIBA-CONICET, Argentina Centre for Mathematical Sciences, Lund University, Sweden
a r t i c l e
i n f o
Article history: Received 25 September 2010 Received in revised form 22 April 2011 Accepted 25 April 2011 Available online 5 May 2011 Keywords: Epidemiology Dengue Individual based model Compartmental model Stochastic
a b s t r a c t We introduce a dengue model (SEIR) where the human individuals are treated on an individual basis (IBM) while the mosquito population, produced by an independent model, is treated by compartments (SEI). We study the spread of epidemics by the sole action of the mosquito. Exponential, deterministic and experimental distributions for the (human) exposed period are considered in two weather scenarios, one corresponding to temperate climate and the other to tropical climate. Virus circulation, final epidemic size and duration of outbreaks are considered showing that the results present little sensitivity to the statistics followed by the exposed period provided the median of the distributions are in coincidence. Only the time between an introduced (imported) case and the appearance of the first symptomatic secondary case is sensitive to this distribution. We finally show that the IBM model introduced is precisely a realization of a compartmental model, and that at least in this case, the choice between compartmental models or IBM is only a matter of convenience. Ó 2011 Elsevier Inc. All rights reserved.
1. Introduction Dengue fever is a vector-born disease produced by a flavivirus of the family flaviviridae [1]. The main vectors of dengue are Aedes aegypti and Aedes albopictus. The research aimed at producing dengue models for public policy use began with Newton and Reiter [2] who introduced the minimal model for dengue in the form of a set of Ordinary Differential Equations (ODE) for the human population disaggregated in susceptible, exposed, infected and recovered compartments. The mosquito population was not modeled in this early work. A different starting point was taken by Focks et al. [3,4] that began by describing mosquito populations in a computer framework named Dynamic Table Model where later the human population (as well as the disease) was introduced [5]. Newton and Reiter’s model (NR) favours economy of resources and mathematical accessibility, in contrast, Focks’ model emphasize realism. These models represent in dengue two contrasting compromises in the standard trade-off in modeling. A third starting point has been recently added. Otero et al. (OSS) developed a dengue model [6] which includes the evolution of the mosquito population [7,8] and is spatially explicit. This last model is somewhat in between Focks’ and NR as it is formulated as a statedependent Poisson model with exponentially distributed times (a detailed description of our model is given in Section 2).
⇑ Corresponding author. Tel.: +46 46 222 0919; fax: +46 46 222 4010. E-mail address: [email protected] (M.A. Natiello). 0025-5564/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2011.04.006
All modeling approaches have been further developed [9–15]. ODE models have received most of the attention. Some of the works explore: variability of vector population [9], human population [10], the effects of hypothetical vertical transmission of dengue in vectors [11], seasonality [12], age structure [13] as well as incomplete gamma distributions for the incubation and infectious times [16]. Contrasting modeling outcomes with those of real epidemics has shown the need to consider spatial heterogeneity as well [17]. The development of computing technology has made possible to produce Individual Based Models (IBM) for epidemics [18,19]. IBM have been advocated as the most realistic models [19] since their great flexibility allows the modeler to describe disease evolution and human mobility at the individual level. When the results are only to be analysed numerically, IBM are probably the best choice. However, they are frequently presented in a most unfriendly way for mathematicians as they usually lack a formulation (expression in closed formulae) and are – at best – presented as algorithms if not just in words [18]. In contrast, working on the ODE side, it has been possible, for example, to achieve an understanding of the influence of distribution of the infectious period in epidemic modeling [20–22]. IBM have been used to study the time interval between primary and secondary cases [23] which is influenced, in the case of dengue, mainly by the extrinsic (mosquito) and intrinsic (human) incubation period. In this work an IBM model for human population in a dengue epidemic is presented. The model is driven by mosquito populations modeled with spatial heterogeneity with the method introduced in [8] (see Section 2). The IBM model is then used to
88
M. Otero et al. / Mathematical Biosciences 232 (2011) 87–95
examine the actual influence of the distribution of the incubation period comparing the most relevant information produced by dengue models: dependence of the probability of dengue circulation with respect to the mosquito population and the total epidemic size. Exponential, delta (fixed times, deterministic) and experimental [24] distributions are contrasted (Section 3). The infectious period and the extrinsic incubation period is modeled using experimental data and measured transmission rates (human to mosquito) [24]. The IBM model produced is critically discussed. We show that it can be mapped exactly into a stochastic compartmental model of a novel form (see Section 4) thus crossing for the first time the valley separating IBM from compartmental models. This result opens new perspectives which we also discuss in Section 4. Finally, Section 5 presents the conclusions of this work.
2. The model It is currently accepted that the dengue virus does not make any effect to the vector. As such, A. aegypti populations are independent of the presence of the virus. In the present model, mosquito populations are produced by the A. aegypti model [8] with spatial resolution of one block using climatic data tuned to Buenos Aires, a temperate city where dengue circulated in the summer season 2008–2009 [25]. The urban unit of the city is the block (approximately a square of 100 m 100 m). Because of the temperate climate the houses are not open as it is often the case in tropical areas. Mosquitoes usually develop in the center of the block which often presents vegetation and communicates the buildings within the block. The model then assumes that mosquitoes belong to the block and not to the houses and they blood-feed with equal probability in any human resident in the block. A. aegypti is assumed to disperse seeking for places to lay eggs. The mosquito population, number of bites per day, dispersal flights and adult mortality information per block is obtained from the mosquito model [8]. The dynamical evolution of the mosquito population is modeled as a Markov Jump Process. It starts with a list of ‘events’ modifying the populations. Each event occurs randomly according to a Kendall–Feller process [26,27] (approximated here as in [28]) and modifies the subpopulations accordingly (e.g. mature eggs turning to larval state, mosquito bite, oviposition, etc.). Mosquitoes can then perform appentential flights and disperse to lay eggs, in contrast to [5] that considers only one of such urban blocks. The model attempts to use only the minimal number of events that are biologically relevant. The events describe basic biological facts and the values of the associated transition rates and their dependence with e.g., temperature are taken from experimental observation. The stochastic framework is motivated by the intrinsic impredictability of natural processes. An interesting issue in the present model is that it has no adjustable parameters, other than the number of breeding sites (BS) per block. Specifically, the model considers six different populations per block: eggs E, larvae L, pupae P, adult females in their first gonotrophic cycle A1, flyer females F and adult females in subsequent gonotrophic cycles A2. The evolution of the six populations in each block is affected by different possible events: death of eggs, egg development and hatching, death of larvae, pupation, death of pupae, pupal development and adult emergence, completion of gonotrophic cycles, death of young adults A1, oviposition by F flyers, death of A2 adults, death of flyers F and flyers dispersal. We consider that during the gonotrophic cycles the mosquito dispersal is negligible, adults A1 and A2 remain near the breeding sites where they grew, until ready to blood feed. After the blood ingest and the completion of their gonotrophic cycles, adult females begin
to fly (becoming flyers F) in search of oviposition sites. Therefore, every time an adult A1 or an adult A2 becomes a flyer one bite was performed. The rates associated to the Markov jump process connecting the different populations are taken from [7,8]. The model disregards differences in virus strains [1], as well as the clinical complications that arise more often in reincident crossinfections. While multiple virus circulation is possible [29,30] its dynamical relevance has yet to be examinated. In this work we attempt to lay the basis for a flexible model that will later develop further incorporating important dynamical effects such as human’s mobility [31], co-circulation of several viral strains and support the testing of proposed public policies, particularly in view that a vaccine for dengue is not expected in the near future [1]. The model calculates every day the mosquito populations and all the events which occur at rates that depend not only on population values but also on temperature, which in turn is a function of time since it changes over the course of the year. Hence, the dependence on the temperature introduces a time dependence in the event rates. The time-step of the model has been fixed at one day. The human population of each block is fixed in the present work and the disease is spatially spread by the mosquito alone. The evolution of the disease in one individual human, h, proceeds as follows: Day d = d0 The virus is transmitted by the bite of an infected mosquito. The transmission probability is set to pmh = 0.75 [2]. Day d = d0 + sE(h) The human h becomes infective (h is said to be exposed to the virus during this period of time). Day d = d0 + sE(h) + j For 1 6 j 6 sI the human h is infective and transmits the virus to a biting mosquito with a probability phm(j). sI indicates the duration in days of the viremic window, in this article taken to be of 5 days following experimental results in [24]; the values of phm(j) are given in Table 1. Day d > d0 + sE(h) + sI The human h is recovered and no longer transmits dengue. The cycle in the human being is then of the form susceptible, exposed, infected, recovered (SEIR). The virus enters the mosquito when it bites a viremic human with a probability phm(j) depending of the day j in the infectious cycle of the specific human bitten (see Table 1). The cycle continues with the reproduction of the virus within the mosquito (extrinsic period), lasting sm days (in this work sm was set to 8 days). After this reproduction period the mosquito becomes infectious and transmits the virus with a probability pmh when it bites. The mosquito follows a cycle susceptible, exposed, infected (SEI) and does not recover [1,5,24,6]. The adult female mosquito population as produced by the A. aegypti simulation is then split into susceptible, sm stages of exposed and one infective compartment according to their interaction with the viremic human population and the number of days elapsed since acquiring the virus. The epidemic starts when one or more humans become viremic. The algorithm followed is: 1. Give individual attributes, sE(h) according to prescribed distribution. 2. Read geometry, size of viremic window, probabilities phm(j), j = 1, . . . , sI, human population in each block, day of the year when the epidemic starts. Table 1 Human-to-mosquito dengue virus transmission probabilities for the different days of the infected stage.
Probability
Day 1
Day 2
Day 3
Day 4
Day 5
0.25
0.8
1
0.95
0.5
89
M. Otero et al. / Mathematical Biosciences 232 (2011) 87–95
6. Repeat all over again from (5a). Each iteration is a new day of the simulation. We take the opportunity to observe that the differences between [6] and the present model are more in the algorithmic side than in the biological content. The present IBM code is substantially more efficient as well as flexible than the compartmental model in [6], allowing for larger cities and better statistics. On
the biological side, the present model is more realistic as it incorporates the experimental values for phm (exponentially distributed in [6]) as well as a more realistic distribution of the extrinsic incubation time.
0.5
Relative frequency
0.4
0.3
0.2
0.1
0
0
1
2
3 4 5 Intrinsic Incubation Period / days
6
7
8
0.4
P
0.6
0.8
1.0
Fig. 1. Distribution for the duration of the exposed period, according to Nishiura’s data [24].
0.0
0.2
N E1 E2 D
0
200
400
600
800
1000
800
1000
0.4
P
0.6
0.8
1.0
BS
0.2
N E1 E2 D
0.0
3. Initialise the blocks with the human population. 4. Read total adult female mosquito population of the day (M), bites, flights to neighbouring blocks and mosquito death probability. Initialise all mosquitoes as susceptible (MS). Set infective bites to zero. 5. Day-loop begins: (a) Calculate the amount of surviving infected mosquitoes. Compute surviving mosquito population with d exposed days, age them by one day. Exposed mosquitoes evolve to infected ones after d = sm days (Use binomial random number generator). (b) Compute probability for a mosquito bite to transmit dengue as pminf = pmhMI/M (compound probability of being infective and being effective in the transmission). Each bite is an independent event according to the underlying mosquito model. (c) Compute the probability for a bite to be made by a susceptible mosquito out of all the non-infective bites potras = MS/ (M(1 pminf)). The non-infective bites have probability (1 pminf), and come from infected mosquitoes failing to transmit the virus, exposed and susceptible mosquitoes, M (1 pminf) = M pmhMI = MS + ME + (1 pmh)MI. (d) Calculate the number of infected and susceptible mosquito bites using binomials with the previous probabilities and the number of total bites. (e) Compute number of humans bitten in each day of the infected state, HI(j). (f) Compute the number of new exposed mosquitoes, taking into account that the probability of human-mosquito contagion is dependent on the stage-day of the human infection. The amount of new infected insects is chosen using binomials. For this purpose, we add the results of the calculation of the number of susceptible mosquitoes that bite humans and PI get infected, M NE ¼ sj¼1 Binðphm ðjÞ; HI ðjÞÞ. Where MNE are the new exposed mosquitoes, Bin(HI(j), phm(j)) is a binomial realization with the day-dependent probability phm(j) and HI(j) the quantity of infected humans bitten by susceptible mosquitoes on infection day j. (g) Perform a random equi-distributed selection of humans bitten by infectious mosquitoes. Build a table of bitten individuals. (h) Update the state of all the humans. If the human belongs to the Susceptible state and has been bitten according to the table built in (5 g) then change the state of selected human to Exposed, record d0 for each exposed human. Susceptible individuals not bitten by an infective mosquito will remain as such and consequently their intrinsic time d remains in 0. Increase the intrinsic time of Exposed and Infected humans by one. Those Exposed individuals for which their intrinsic time has surpassed the value d = d0 + sE(h) are moved to the Infected state, while those Infected individuals whose intrinsic time is larger than d = d0 + sE(h) + sI, are moved to the Removed state (i) Compute the number of individuals in each state for every cell. (ii) Read the total adult female mosquito population of the day (M), bites, flights to neighboring blocks.
0
200
400
600
BS Fig. 2. Probability of local circulation of virus (probability of having at least one secondary case). Top: with a tropical temperature scenario, Bottom: with a temperate climate, BS: Number of breeding sites per block.
90
M. Otero et al. / Mathematical Biosciences 232 (2011) 87–95
3. Epidemic dependence on the distribution of the exposed period We implemented four different distributions for the duration of the exposed period assigned to human individuals: Nishiura’s experimental distribution [24], a delta and exponential distributions with the same mean that the experimental one and an exponential distribution with the same median than the experimental one. We call them N, D, E1 and E2, respectively. See Fig. 1 for an histogram of Nishiura’s experimental distribution. The study was performed in two different climatic scenarios, one with constant temperature of 23 degrees Celsius, that represents tropical regions and one with the mean and amplitude characteristic of Buenos Aires, a city with temperate climate. The number of effective breeding sites [7] was varied between 50 and 1000. The simulations were performed using identical mosquito populations in all cases (same data file), for an homogeneous urban area of 20 by 20 blocks, hosting 100 people per block, making a total of 40,000 human beings. Hence, all differences correspond to the disease dynamics that was previously identified as the main source of stochastic variations. The simulations with temperate climate were started on January 1st, i.e., ten days after the summer solstice (December 21st in the southern hemisphere). The reported statistics is computed by averaging the outcomes of 1000 simulations with different seeds for the pseudo-random routines. The mosquito part of the simulation has been successfully tested in larger modeling environments, suggesting that just mosquito dispersal is an important ingredient contributing to the mosquito
persistence (survival) in cities with temperate climate [8]. The human population in this model is taken to be static, with fixed population density, hence a larger population size will not change the overall results. The effects of moderately changing the human population can be easily understood. The window for secondary cases will not change (does not depend on the population), the epidemic size will be proportional to the human population size, while the probability of virus-circulation as well as R0 decrease when the population increases (the probability of the index case to be bitten by a mosquito varies as the inverse of the population size [32]), we therefore do not present confirmatory results of these well-understood effects. The main questions were: considering the total number of recovered individuals, how is the distribution of epidemic sizes influenced by the choice of distribution? How does the probability for having no secondary cases change? How is the predicted duration of the epidemic influenced by our choices? And finally, how is the distribution of time between epidemiologically related cases affected? This time is called the generation time in [33]. The epidemiologically relevant question is: can a given (dengue) case be epidemiologically related to (one of) the previously detected case(s)? The proper answer to the later question will be of the type: considering the time interval separating both events, the probability of they being epidemiologically related is . . . Generation times need to be studied at the probability distribution level. We will present the minimal information of the distributions in the form of box-plots, showing the median, lower and upper quartiles, and the ‘outlier points’. Before showing our results, it is worth to realise that the probability of having no secondary cases will not be sensitive
E1
30000
30000
30000
0 00
250
0 200
10000
150
10000
100
20000
50
20000
10
Epidemics Size
30000
50
250
E2 40000
10
250
D
Epidemics Size
200
BS
40000
BS
200
BS
150
10
150
250
0 200
0 100
10000
50
10000
150
20000
100
20000
50
Epidemics Size
40000
10
Epidemics Size
N 40000
BS
Fig. 3. Box plot graphs for the epidemic size (total number of infected humans) at constant temperature. Top-left: N-distribution, top-right: E1-distribution, bottom-left Ddistribution and bottom-right: E2-distribution. BS: Number of breeding sites per block.
91
M. Otero et al. / Mathematical Biosciences 232 (2011) 87–95
to the choice of distribution for the exposed period, as this probability depends only on the probability of the introduced case being bitten by the mosquitoes, the probability of the mosquitoes of acquiring the virus, surviving the extrinsic period and finally transmitting the virus to a human in a bite. Nothing in this process depends on the choice of distribution. In contrast, we expect the distribution of times between epidemiologically related cases to depend strongly on the choice of distribution, since it reflects the sum of the two incubation periods (intrinsic and extrinsic). Some intuition of the epidemiological situation described can be achieved calculating R0, the basic reproductive number. We follow the definition in [32] and define R0 as the average number of secondary cases produced by a single case when the epidemic starts. A rapidly growing epidemy will have R0 > 1 while a disease going to extinction will display R0 < 1. A plot of this R0 as a function of the number of breeding sites is displayed in Fig. 5. In the two conditions studied in this manuscript (summer conditions in temperate climate and a year-long favourable tropical climate), R0 > 1 at values of the breeding sites per block as low as 15. Confirming the reasoning above, there is no sensitivity for the probability of having local circulation of the virus (defined as having at least one secondary case following the introduced case) with the statistical distribution of the exposed period, see Fig. 2. However, we note that our model does not force local circulation of the virus, since the probability of circulation is not identically equal to one. According to the Law of Large Numbers, in a sufficiently large number of simulations, a portion of them will yield extinction of the disease.
The size of the epidemic at constant temperatures makes a transition from very small outbreaks to a large outbreak reaching almost all the people in the simulation. The transition happens in the region 50–100 breading sites per block for the four distributions studied, see Fig. 3. We detected no epidemiologically important differences produced by the use of one or another distribution function. The size of epidemics with seasonal dependence is presented in Fig. 4. The results for the D-distributions and N-distribution do not present differences. The E1-distribution (equal mean) overestimates the final size while the E2-distribution substantially agrees with the N (experimental) one. Both exponential distributions present a larger variance than the N-distribution. This difference may matter when worst-possible scenarios are considered. The size of the epidemic outbreak begins to increase with the number of breeding sites in the 100–150 region for all the distributions. It is possible to argue that the differences in epidemic size observed for different distributions in the case with seasonal dependence correspond to a faster evolution of the epidemic outbreak during the time-window of favourable conditions. On the other hand, in the constant temperature scenario the epidemic outbreak stops because of the decrease in susceptible people produced by the epidemic. In this case we observe no significant differences, see Fig. 6, the epidemic size always ends up around 40,000 individuals. However, the epidemics evolves faster the larger the number of breeding sites (rightmost plots). It is worth observing that major outbreaks last more than one year in this 40,000 people urbanization. Since there is no human movement incorporated, the dura-
25000
BS
900
950
1000
950
1000
850 850
900
750
800 800
700 700
750
600
650
600
650
500
550
500
450
400
350
300
250
200
950
1000
900
850
800
750
700
650
600
550
500
450
0 400
0 350
5000
300
5000
150
10000
100
10000
15000
50
15000
20000
10
Epidemics size
20000
550
E2
25000
250
450
D
25000
200
350
BS
30000
150
400
BS
30000
100
250
10
950
1000
850
900
800
750
650
700
550
600
500
400
450
300
350
200
250
150
0 100
0 10
5000
50
5000
300
10000
200
10000
15000
150
15000
20000
50
20000
100
Epidemics size
25000
10
Epidemics size
E1 30000
50
Epidemics size
N 30000
BS
Fig. 4. Box plot graphs for the epidemic size with seasonal dependence. Top-left, N-distribution, top-right, E1-distribution, bottom-left D-distribution and bottom-right, E2distribution. BS: Number of breeding sites per block.
92
M. Otero et al. / Mathematical Biosciences 232 (2011) 87–95
30
tious material, and the appearance of symptoms. . . Several implications can be drawn from these epidemiological principles and serve as a conceptual model for an individual infection process. First, it is most appropriate to represent the individual infection as a series of discrete events and periods. Second, the discrete events have no duration in their own right and only trigger the change between infection periods. Third, the discrete periods indicate the infection status of an individual and are part of the individual’s characteristics.
D N E1 E2
25
R0
20
15
10
5
0
0
200
400
600
800
1000
800
1000
BS
30
D N E1 E2
25
R0
20
15
10
5
0
0
200
400
600
BS Fig. 5. Estimation of the basic reproductive number R0 as a function of the number of breeding sites per block. Top: with a tropical temperature scenario. Bottom: with a temperate climate.
tion of the outbreak depends critically on the dispersion of female mosquitoes. At this point we could ask: is there any relevant statistics that depends on the distribution of exposed times? The answer is yes. Assume a case of imported dengue is detected, how long do we have to wait to know if we are facing an outbreak or not? The time elapsed between the primary and the first secondary case is given by adding the extrinsic incubation period and the exposed time. Box plots produced after 1000 simulations of outbreaks wit a mosquito population supported on 200 BS per block are displayed in Fig. 7. We observe that in this case the exponential distributions exaggerate the dispersion of results producing too early as well as too late cases as compared with the experimental distribution while the delta-distribution compresses the ‘alert window’ too much with respect to the experimental distribution. 4. Is the IBM model an implementation of a compartmental model? The use of IBM in epidemiology is usually advocated on an ontological perspective [34,19], we quote the argument in [19] . . .the epidemiology literature has always described an infection history as a sequence of distinct periods, each of which begins and ends with a discrete event. The critical periods include the latent, the infectious, and the incubation periods. The critical events include the receipt of infection, the emission of infec-
The evolution of dengue at the individual level is described in detail in the literature, including experimental results [24] which are seldom available for other illnesses. Thus, dengue is a good case to put the thesis at test. The description in terms of events immediately calls for stochastic population models, while the different human to mosquito transmission probabilities can be handled easily introducing age structure in the infective human population, The description of the proper probability distribution for the exposed (latent) period is the major obstacle towards a compartmental stochastic model. Each human individual spends k days, k 2 {1, . . . , sE}, in the exposed state before becoming infective. The probability of any individual to spend k days is P(k). Hence, in the group of N individuals that become infected at d = d0, a portion N(k) will spend k days before they become infective, where N(k) is a random deviate taken from Multinomial(N, P(1), . . . , P(sE)), a multinomial distribution. This is the only information available to the simulation. IBM generate additional structure, since each individual is assigned to one of the classes N(1), . . . , N(sE) in an entirely random way. This apparent additional information is in fact arbitrary because of the randomness and it is averaged out when presenting the results. Although the algorithm of the IBM model assigns the time spent in the exposed class to the individual, the final cause of having a distribution of exposed times is not known (neither to us nor to the algorithm). Different exposed times may arise either because of differences (in virus resistance) in the exposed individuals or because of differences in the incoming virus due e.g., to biological processes concerning the development and transmission of the virus by the mosquito. Let us define
b PðkÞ ¼
sE X
PðjÞ;
ð1Þ
j¼k
QðkÞ ¼
PðkÞ ¼ Pðj ¼ k=j P kÞ; b PðkÞ
bk ¼ N N
k X
NðjÞ:
ð2Þ
ð3Þ
j¼1
Clearly Q(k) 6 1 and Q(sE) = 1. Consider the numbers E(1) = N, E(k) = Binomial(E(k 1), 1 Q(k)) and N(k) = E(k 1) E(k) for 2 6 k 6 sE. We will now recall a few elementary results regarding multinomial distributions. We write Multi for the multinomial distribution. Lemma 1. Multinomial splitting
MultiðN; pð1Þ; . . . ; pðk 1Þ; pðkÞ; . . . ; pðsE ÞÞ ¼ b ¼ Multi N; pð1Þ; . . . ; pðk 1Þ; PðkÞ b Multi NðkÞ; Q ðkÞ; . . . ; Q ðsE Þ :
ð4Þ
Proof. The proof is simple algebra simplifying the expression for the probabilities in the right side of the equation. h
93
M. Otero et al. / Mathematical Biosciences 232 (2011) 87–95
0 500
1500
2500
3500
0
200
400
600
800
0
200
400
600
800
0
200
400
600
800
0 500
1500
2500
3500
0
200
400
600
800
0
200
400
600
800
0
200
400
600
800
0 500
1500
2500
3500
0
200
400
600
800
0
200
400
600
800
0
200
400
600
800
0 500
1500
2500
3500
0
200
400
600
800
0
200
400
600
800
0
200
400
600
800
Fig. 6. Histogram of the duration (in days) of the epidemic outbreak at constant temperature for different number of breeding sites, left to right 50, 100, 150 and 200 breeding sites (BS) per block. Top line, N-statistics, second line E1-statistics, third line E2-statistics and bottom line D-statistics. x-axis: Duration (in days); y-axis: Relative frequency.
Corollary 1 (Binomial decomposition). The multinomial distribution can be fully decomposed in terms of binomial distributions in the form
MultiðN; pð1Þ; . . . ; pðsE ÞÞ ¼
sE Y
BinomialðEðkÞ; QðkÞÞ:
ð5Þ
k¼1
Proof. Repeated application of the previous lemma is all what is needed. h We now state the application of these results to our modeling problem. Theorem 1. Compartmental presentation Consider sE compartments associated to the days k = 1, . . . , sE after receiving the virus from the mosquito for those humans that are not yet infective. Let the population number in the k compartment be E(k), and the number of humans that become infective on day k be N(k) which is a random deviate distributed with Binomial(E(k), Q(k)), with the definitions given above. Then the exposed period that corresponds to the individuals is distributed with P(k).
Proof. It follows immediately from the corollary.
h
The resulting scheme for the compartmental model is illustrated in Fig. 8. 5. Summary, discussion and conclusions We have developed an IBM model for the evolution of dengue outbreaks that takes information from mosquito populations sim-
ulated with an A. aegypti model and builds thereafter the epidemic part of the evolution. The model was used to explore the actual influence of the distribution of exposed time for humans in those characteristics of epidemic outbreaks that matter the most: determining the level of mosquito abundance that makes unlikely the occurrence of a dengue outbreak and determining the size and time-lapse of the outbreak. The distributions used are (a) an experimentally obtained distribution (Nishiura [24]) (labeled N), (b) and (c) exponential distributions adjusted to have the same mean, respectively the same median as N, labeled E1 and E2, and (d) a fixed-time distribution equal to the experimental mean (labeled the D-distribution). 5.1. Main findings For mosquito populations sufficiently large to support the epidemic spread of dengue, most of the stochastic variability is provided by the epidemic process rather than by stochastic fluctuations in the mosquito population. Thus, the splitting of the models improves substantially the performance of the codes. The probability of producing one or more secondary cases after the arrival of an infective human does not depend on the choice of distribution. The characteristic size of the epidemics under a temperate climate are exaggerated by the E1-distribution but presents no substantial difference for the other distributions. The dispersion of values is exaggerated by both exponential distributions. We observed no important differences in the duration of epidemics developed under a constant temperature since the outbreaks reach almost all the population and the velocity is regulated by the dispersion of the mosquitoes in the absence of movement by humans.
N
E1
60 30 20 10 0
0
0
10
10
20
20
30
30
40
40
50
50
60 50 40
50 40 30 0
10
20
time / days
60
M. Otero et al. / Mathematical Biosciences 232 (2011) 87–95
60
94
E2
D
Fig. 7. Distribution of times between the arrival of an infectious person and the time for the first secondary case. The scenario corresponds to 200 breeding sites per block at constant temperature of 23 °C and 100 people per block. Results correspond to the delta-distribution (D), the experimental distribution (N), the exponential E2 and E1 distributions.
S
E(1)
1−Q(1)
Q(1)
E(2)
E(Te)
Q(Te)=1
Q(2)
Pmh I(1)
1
I(2)
I(Ti)
1 R
Phm(1) 1 IM
EM(Tm)
EM(1)
SM
Phm(Ti)
Fig. 8. Scheme for the progression of dengue. Circles for human subpopulations (susceptible, exposed by day, infective by day and recovered) and squares for mosquito subpopulations. Solid directed arrows indicate daily progression, dotted directed arrows indicate several days of progression, bi-directed arrows with a circle indicate interactions resulting in transitions for members of one population. Probabilities are indicated next to arrows. The mosquito dynamics (birth, death, . . . ) is not represented.
The only statistic able to discriminate easily between the four distributions of exposed time was found to be the time of appearance of the first secondary case (generation time). A result that was expected as well. This fact allows to conjecture that in an extreme situation for the mosquito development, e.g., approaching winter temperatures in a temperate climate, the different delays given by the four distribution will discriminate different behaviours. In the limit, if the incubation time for the first infected case were so long that there are no mosquitos around to pick up the virus, then no disease circulation will occur. That could only happen for extremely late (seasonally speaking) primary infections. 5.2. Limitations The split between mosquito evolution and epidemic evolution is not perfect since the events bite and flight are treated as independent events while they are in fact correlated [35–38], both being related to oviposition.
A general feature of stochastic models is the fact that populations are subdivided into compartments. Within each class, all individuals behave as ‘average individuals’. Extensive model-testing and falsation attempts are the only tool to gain confidence in the proposed structure.
5.3. Concluding remarks The model we developed is shown to be an IBM implementation of a compartmental model, of a form not usually considered. At the level of the description in this work, IBM does not play a fundamental role, contrary to what it has been previously argued [34,19]. Rather, its use is a matter of algorithmic convenience. This result indicates that ‘IBM versus compartmental models’ is not a fundamental dichotomy but it may be a matter of choice (depending of the skills and goals of the user): IBM facilitate coding, compartmental models lend themselves to richer forms of analysis.
M. Otero et al. / Mathematical Biosciences 232 (2011) 87–95
In conclusion, only very specific matters seem to depend on the characteristics of the distribution of exposed times for human beings. The question: do the present results support the inclusion of highly detailed incubation periods? cannot be answered without further specification. The statistical detail required in a given simulation depends on the questions that we intend to answer with the model and on the use to be given to these answers. In the present work the answer to the question: what could go wrong if we omit (or misrepresent) the distribution of exposed times? has been explored. For example, an accurate description of the distribution of exposed times is mandatory to assess if a (supposedly secondary) dengue case is epidemiologically related to a previous one. Although other features of the dynamics appear to be less sensitive to this distribution, we cannot expect this always to be the case. Every omission represents a potential liability but on the other hand, in any conceivable model, detail is always finite. Perhaps one of the most striking and unexpected results of the present work is that the exponential distribution with the same median as the experimental data is to be preferred in front of the one with the same mean, contrary to the popular belief (which was indeed what we believed when we started the work). Looking towards the past, conclusions reached using exponential and delta distributions cannot be objected on such basis. Looking towards the future, simple compartmental models can be constructed as well using realistic distributions and there is no reason to limit the models to the choice: exponential, gamma or delta. References [1] D.J. Gubler, Dengue and dengue hemorrhagic fever, Clin. Microbiol. Rev. 11 (1998) 480. [2] E.A.C. Newton, P. Reiter, A model of the transmission of dengue fever with an evaluation of the impact of ultra-low volume (ulv) insecticide applications on dengue epidemics, Am. J. Trop. Med. Hyg. 47 (1992) 709. [3] D.A. Focks, D.C. Haile, E. Daniels, G.A. Moun, Dynamics life table model for Aedes aegypti: Analysis of the literature and model development, J. Med. Entomol. 30 (1993) 1003. [4] D.A. Focks, D.C. Haile, E. Daniels, G.A. Mount, Dynamic life table model for Aedes aegypti: Simulations results, J. Med. Entomol. 30 (1993) 1019. [5] D.A. Focks, D.C. Haile, E. Daniels, D. Keesling, A simulation model of the epidemiology of urban dengue fever: Literature analysis, model development, preliminary validation and samples of simulation results, Am. J. Trop. Med. Hyg. 53 (1995) 489. [6] M. Otero, H.G. Solari, Mathematical model of dengue disease transmission by Aedes aegypti mosquito, Math. Biosci. 223 (2010) 32. [7] M. Otero, H.G. Solari, N. Schweigmann, A stochastic population dynamic model for Aedes aegypti: Formulation and application to a city with temperate climate, Bull. Math. Biol. 68 (2006) 1945. [8] M. Otero, N. Schweigmann, H.G. Solari, A stochastic spatial dynamical model for Aedes aegypti, Bull. Math. Biol. 70 (2008) 1297. [9] L. Esteva, C. Vargas, Analysis of a dengue disease transmission model, Math. Biosci. 150 (1998) 131. [10] L. Esteva, C. Vargas, A model for dengue disease with variable human population, J. Math. Biol. 38 (1999) 220. [11] L. Esteva, C. Vargas, Influence of vertical and mechanical transmission on the dynamics of dengue disease, Math. Biosci. 167 (2000) 51. [12] L.M. Bartley, C.A. Donnelly, G.P. Garnett, The seasonal pattern of dengue in endemic areas: Mathematical models of mechanisms, Trans. R. Soc. Trop. Med. Hyg. 96 (2002) 387. [13] P. Pongsumpun, I.M. Tang, Transmission of dengue hemorrhagic fever in an age structured population, Math. Comput. Model. 37 (2003) 949.
95
[14] K. Magori, M. Legros, M.E. Puente, D.A. Focks, T.W. Scott, A.L. Lloyd, F. Gould, Skeeter buster: A stochastic, spatially explicit modeling tool for studying Aedes aegypti population replacement and population suppression strategies, PLoS Negl. Trop. Dis. 3 (9) (2009) e508, doi:10.1371/journal.pntd.0000508.