Spatio-temporal characteristics of dengue outbreaks
Saulo D. S. Reis, Lucas Böttcher, João P. da C. Nogueira, Geziel S. Sousa, Antonio S. Lima Neto, Hans J. Herrmann, José S. Andrade Jr
SSpatio-temporal characteristics of dengue outbreaks
Saulo D. S. Reis, ∗ Lucas B¨ottcher,
2, 3, 4, † Jo˜ao P. da C. Nogueira, Geziel S.Sousa,
5, 6
Antonio S. Lima Neto,
5, 7
Hans J. Herrmann,
1, 8 and Jos´e S. Andrade Jr. Departamento de F´ısica, Universidade Federaldo Cear´a, 60451-970 Fortaleza, Cear´a, Brazil Computational Medicine, UCLA, 90095-1766, Los Angeles, United States Institute for Theoretical Physics, ETH Zurich, 8093, Zurich, Switzerland Center of Economic Research, ETH Zurich, 8092 Zurich, Switzerland Secretaria Municipal de Sa´ude de Fortaleza (SMS-Fortaleza), Fortaleza, Cear´a, Brazil Departamento de Sa´ude Coletiva, Universidade Federal do Cear´a, Fortaleza, Cear´a, Brazil Universidade de Fortaleza (UNIFOR), Fortaleza, Cear´a, Brazil PMMH. ESPCI, 7 quai St. Bernard, 75005 Paris, France (Dated: June 5, 2020)
Abstract
After their re-emergence in the last decades, dengue fever and other vector-borne diseases area potential threat to the lives of millions of people. Based on a data set of dengue cases in theBrazilian city of Fortaleza, collected from 2011 to 2016, we study the spatio-temporal characteristicsof dengue outbreaks to characterize epidemic and non-epidemic years. First, we identify regionsthat show a high prevalence of dengue cases and mosquito larvae in different years and also analyzetheir corresponding correlations. Our results show that the characteristic correlation length ofthe epidemic is of the order of the system size, suggesting that factors such as citizen mobilitymay play a major role as a drive for spatial spreading of vector-borne diseases. Inspired by thisobservation, we perform a mean-field estimation of the basic reproduction number and find thatour estimated values agree well with the values reported for other regions, pointing towards similarunderlying spreading mechanisms. These findings provide insights into the spreading characteristicsof dengue in densely populated areas and should be of relevance for the design of improved diseasecontainment strategies. ∗ saulo@fisica.ufc.br † [email protected] a r X i v : . [ phy s i c s . s o c - ph ] J un NTRODUCTION
According to a recent report of the world health organization (WHO), over 40% of theworld’s population are at the risk of a dengue infection [1]. Alarmingly, in the last half-century, a 20 to 30-fold increase of dengue cases has been monitored world-wide as shownin Fig. 1 (left). Dengue fever is a vector-borne disease and thus requires a disease vectorto transmit the virus from one infected human to another susceptible one. Dengue virustransmission occurs through bloodsucking female
Aedes aegypti mosquitoes. Infected hu-mans transmit the virus, up to maximally 12 days after first symptoms occur, to a mosquitowhich further transmits the virus after an incubation period of 4-10 days [1]. Symptomsinclude high fever, headache, vomiting, skin rash, and muscle and joint pains [2]. Aedesaegypti is also responsible for the transmission of other severe vector-borne diseases suchas yellow fever, chikungunya and zika whose recent outbreaks are challenging health offi-cials in different countries [3–5]. Vector-borne diseases are responsible for the death of morethan one million humans each year, disrupt health systems, and obstruct the developmentof many countries [1]. As a consequence of the unavailability of vaccines and/or drug re-sistance, as is the case for many vector-borne diseases, control measures such as personalprotection, reduction of vector breading habitats, and usage of insecticides are essential tocontain outbreaks [1, 6].Recently, studies have demonstrated the poor correspondence between levels of infesta-tion measured by entomological surveys at national and local level [8–10]. The correlationbetween high infestation, measured by larval entomological indicators, and risk of dengueepidemic has been shown to be weak in the most diverse scenarios [11–15]. Despite theevidence of the low usefulness and incipient positive predictive value of these larval surveysbeing extensively demonstrated, some countries still consider the thresholds of infestationindices as a main indicative of epidemic risk and trigger of spatially-oriented vector controlactions [8, 9, 15].In fact, there is pessimism about the ability of larval indicators to be good parametersfor the prediction of epidemics. Several authors have suggested the systematic use of pupalindices and those that directly measure the density of the adult mosquito, even if thisimplies a radical reformulation in the routines of vector control programs. They argue thatthe larval indexes are no longer able to fulfill their main objective, which is, when estimated2 igure 1.
Temporal evolution of the number of dengue cases.
In the left panels, we showthe number of dengue cases that were reported to the WHO during the last half-century [1, 7].Between 1995–2007, the shown numbers are averages over the indicated periods. The number ofcases has been increasing drastically during the last half-century. In the right panel, we show thenumber of reported dengue cases per day in Fortaleza from 2011 to 2016. the abundance of the adult vector, to anticipate the transmission of the disease [8, 10, 13,15–17]. Moreover, pupal surveys and mosquito capture are of low cost-effectiveness andinaccurate in low infestation scenarios where the mosquito has high survival due to optimalclimatic conditions and other factors, increasing vector competence and causing sustainedtransmission of the virus with low pupal and mosquito densities.Urban mobility, as expressed by the movement of viremic humans, has become a more im-portant factor than levels of infestation and abundance of the Aedes aegypti vector [18, 19].Several studies point out the importance of the human movement in the spatio-temporaldynamics of urban arboviruses transmitted by Aedes aegypti, in particular, dengue fever [20–22]. In 2019, the Pan American Health Organization (PAHO) published a technical doc-ument that points out that the risk stratification of transmission, using the most diverseinformation available, is the path to better efficiency of vector control [23]. Theoretically,mapping cities in micro-areas with different levels of transmission risk could make vectorcontrol strategies more effective, especially in regions with limited human and financialresources [24–26]. However, human mobility patterns are still not subject to accurate moni-toring and can be crucial for the definition of high-risk spatial units that should be prioritizedby the national and local programs for the control of urban arboviruses transmitted by Aedes3egypti [8].Predicting and containing dengue outbreaks is a demanding task due to the complexnature of the underlying spreading process [27, 28]. Many factors such as environmentalconditions or the movements of humans within certain regions influence the spreading ofdengue [29–32]. The increasing danger of dengue infections makes it, however, necessaryto foster a better understanding of the spreading dynamics of dengue. In this work, westudy the spatio-temporal characteristics of dengue outbreaks in Fortaleza from 2011-2016to characterize epidemic and non-epidemic years. Fortaleza, the capital of Cear´a state isone of the largest cities in Brazil and is located in the north-east of the country wheredengue and other neglected tropical diseases show a high prevalence [16, 33]. We showin Fig. 1 (right) that up to 1000 dengue cases have been reported per day in Fortalezaduring 2011 and 2016. In our analysis, we identify regions that exhibit a large number ofdengue infections and mosquito larvae in different years, and also analyze the correspondingcorrelations throughout all neighborhoods of Fortaleza. We observe that the correlationlength (i.e., the characteristic length scale of correlations between case numbers at differentlocations in the system) is of the order of the system size. This provides strong support tothe fact that presence factors such as citizen mobility driving the spatial spreading of thedisease.Motivated by the observation that people interact across long distances, we use a mean-field model to compare the outbreak dynamics of two characteristic epidemic years withcorresponding analyses made for other regions [34]. In particular, we perform a BayesianMarkov chain Monte Carlo parameter estimation for a mean-field susceptible-infected-recovered (SIR) model which has been previously successfully applied to the modelingof dengue outbreaks [34]. Our results provide insights into the spatio-temporal character-istics of dengue outbreaks in densely populated areas and should therefore be relevant formaking dengue containment strategies more effective.
DATA SET
In total, we consider three data sets that are available as
Supplemental Material . The firstdata set consists of spatio-temporal information on dengue infections in the city of Fortalezafrom 2011 to 2016. These data have been provided by the
Epidemiological Surveillance ear Number of cases reported Number of dengue cases reported in Fortaleza between 2011 and 2016.
Division of Fortaleza Health Secretariat , and contains both the date when an infected per-son reports a potential dengue infection to a physician and the corresponding geographiclocation. Serology and clinical diagnosis were used to diagnose dengue infections. Datasetthat was provided by Fortaleza Health Secretariat has been anonymized. Epidemiologicaland clinical variables have also been removed. We remove cases of infected tourists andpeople who live in the metropolitan area of Fortaleza. The total number of dengue casesbetween 2011 and 2016 in our data set after data processing is 133540. Table I shows thenumber of cases reported for each year. Our second data set contains information aboutAedes aegypti larvae measurements in Fortaleza from 2011 to 2016. Measurement data isavailable in fortnight intervals (i.e., every two weeks) from 113 strategic points (SP) (e.g.junk yards) that are monitored by health authorities of Fortaleza over the whole city. ASP is said to be positive if Aedes aegypti larvae have been found independent of the actualamount. Data on Aedes aegypti infestation according to each Strategic Point are availableon the Fortaleza Daily Disease Monitoring System [35] and were tabulated and consolidatedby the epidemiological surveillance division team. Due to the influence of precipitation onthe mosquito population size [29], we also consider data from three different rain gauges inFortaleza during 2011 to 2016.
Dengue outbreaks in Fortaleza
According to the Brazilian ministry of health and the health administration of the prefec-ture of Fortaleza, the population of Fortaleza had been at a high risk of a dengue infection5
Cycle M ea s u r e Dengue casesRainfall ( × ×
10) 0 4 8 12 16 20 24
Cycle M ea s u r e Dengue casesRainfall ( × × Cycle M ea s u r e Dengue casesRainfall ( × ×
10) 0 4 8 12 16 20 24
Cycle M ea s u r e Dengue casesRainfall ( × × Cycle M ea s u r e Dengue casesRainfall ( × ×
10) 0 4 8 12 16 20 24
Cycle M ea s u r e Dengue casesRainfall ( × Figure 2.
The number of reported dengue cases, positive SPs and rainfall as a functionof time from 2011 to 2016.
For each year from 2011 until 2016, we show the number of reporteddengue cases, positive SPs and the rainfall as a function of time. One time interval (cycle) is onefortnight. The number of positive SPs and amount of rain in liters per square meter have beenrescaled by a factor of 10. There is no larvae measurement data available for 2016. during the years 2011, 2012, 2015, and 2016 that are classified as epidemic years [36]. Thenumber of reported dengue cases in this time period is given in Table I and shown in Fig. 1(right). Considering the population size of 2.5 million [37], an alarmingly high number ofseveral hundred up to almost 1000 new infections per day have been reported during thistime period. As shown in Fig. 2, the number of reported dengue cases starts to increase6 orrelations Year τ max (D , R) 5 3 2 2 2 0 τ max (D , SP) 2 2 1 3 3 - τ max (SP , R) 2 1 -1 -1 -1 -Table II.
The maximum correlation coefficient time lags from 2011 until 2016.
Foreach year from 2011 until 2016, we compute the correlation coefficients of the number of reporteddengue cases (D), positive SPs and the rainfall (R) for different time lags τ . The value of τ = 1corresponds to one fortnight. We show the time lags that correspond to a maximum correlationcoefficient. There is no SP data available for 2016. in January and February just shortly after the beginning of the rain season and typicallyreaches its peak before July. This shows that the corresponding climate conditions facilitatethe growth of mosquito populations. In addition, we also show the time evolution of thenumber of positive SPs (i.e., the number of strategic measurement points where Aedesaegypti larvae have been found). The data in Fig. 2 indicates that the number of positiveSPs changes according to the rainfall, because the curves behave very similar although theirrelative amplitudes vary substantially from year to year.To better understand the correlations between rainfall, larvae growth, and dengue caseswe compute the corresponding correlation coefficients for different time lags τ . A time lagof τ = 1 means that the two considered curves are shifted by one fortnight. In Table II,we show the time lags τ max that correspond to the maximum correlation coefficient. Inthe case of dengue occurrences and rainfall, we find a mean value of ¯ τ max ( D, R ) = 2 . τ max ( D, SP ) = 2 . igure 3. Dengue outbreaks in Fortaleza from 2011 until 2016.
Heat maps of all denguecases in Fortaleza from 2011 until 2016. Blue areas correspond to regions with a low prevalenceof dengue cases whereas red ones indicate large dengue outbreaks. The epidemic years are 2011,2012, 2015 and 2016. All heat maps have been generated with Folium [38]. number of positive SPs is found after the rainfall maximum, whereas the opposite situationoccurs in other years, which means that they are uncorrelated. This can be also seen bycomparing the time evolutions of rainfall and the number of positive SPs in Fig. 2. Apossible reason to explain this result may be the fact that a SP was considered as positiveeven if a small amount of larvae was found at this location. For this reason the positivenessof a SP indicates the spatial extent of mosquitoes and not their concentration. Therefore,correlations between rainfall and positive SPs have to be interpreted with caution. Still,the very similar shapes of the positive SP and the rainfall curves indicate that the spatialspread of the mosquito changes with the amount of rain. The influence of climate conditionson dengue outbreaks in Brazil has also been recently analyzed in reference [29]. Due tothe equatorial location of Fortaleza, the mean temperature over one year is 26 . ± . ◦ Cand can be regarded as constant with only limited impact on the local dengue spreading8 igure 4.
Positive SPs in Fortaleza from 2011 until 2015.
Heat maps of the larvae measure-ment outcomes in Fortaleza from 2011 until 2015. Each district has a different number of SPs. Iflarvae are found, the SP measurement is said to be positive. Blue areas correspond to regions witha low number of positive SPs, whereas red ones indicate the opposite. All heat maps have beengenerated with Folium [38]. dynamics [39].In addition to the temporal information on dengue cases in Fortaleza, we collected dataabout the geographical locations of all infection cases. We illustrate the geographical dis-tribution of dengue cases in Fig. 3 from 2011 until 2016. If the total number of dengueinfections in a certain area is large over the course of one year, this area appears in red.Areas with small numbers of dengue infections are colored blue. In 2011 and 2012 the mainoutbreaks occur in similar regions and cover almost the whole city except some parts in theeastern outskirts. In the non-epidemic years 2013 and 2014, the overall dengue prevalence isclearly reduced. However, some parts in the city center in the north west that also exhibitlarge numbers of dengue cases in the epidemic years are still at a high risk of dengue out-breaks. Also in the epidemic years 2015 and 2016, some neighborhoods in the south show a9igh dengue prevalence.In addition to the geographical location of dengue cases, we also show a heat map of thespatial distribution of positive SPs in Fig. 4. Interestingly, and in contrast to the spatialdistribution of dengue cases, the majority of positive SPs are located at the same place inFortaleza regardless of the year. In particular the city center and some parts in the northand south west of the town exhibit a large number of positive SPs. While according toFig. 4 the regions of major outbreaks change substantially between 2011, 2013 and 2015,the number of positive SPs does not at all. Concerning the outbreak years 2015 and 2016,the recurrent dengue outbreaks and numbers of positive SPs seem to be even anti-correlated.These observations suggest that an effective disease intervention measure should target theneighborhoods which exhibit recurrent dengue outbreaks instead of numbers of positive SPs.In particular, the density of SPs is not homogeneous, which could be another reason whythey fail to reproduce the major outbreaks of 2015 and 2016 in the south east of the city.
SPATIO-TEMPORAL CORRELATIONS
For a better understanding of the spreading dynamics of dengue in urban environments,it is crucial to analyze the influence of local dengue outbreaks on outbreaks at other locationsand vice-versa. In this way, We therefore now focus on spatio-temporal correlations betweendengue incidences of different districts. In total, there are 118 neighborhoods in Fortaleza,and for each neighborhood i we consider the number of reported dengue cases NC ti at time t and the corresponding population POP i . The disease prevalence in neighborhood i at time t is given by p ti = NC ti POP i . (1)We characterize the spatio-temporal correlations of dengue outbreaks between sites i and j by c ij ( τ, r ij ) = 1 T (cid:80) Tt =1 (cid:0) p t − τi − (cid:104) p i (cid:105) (cid:1) (cid:0) p tj − (cid:104) p j (cid:105) (cid:1) σ σ , (2)where τ denotes a time lag, r ij is the distance between neighborhood i and neighborhood j , and σ i = T − (cid:80) Tt =1 ( p ti − (cid:104) p i (cid:105) ) represents the variance. In our following analysis, weconsider two time intervals: (i) four weeks (i.e., T = 12) and (ii) two weeks (i.e., T = 26).Dengue outbreaks that occur in a certain region of the city may lead to outbreaks in10 r i j . . . . . . h c ij ( τ ) i − . − . . . . c i j ( τ ) P ( c ij ) r i j . . . . . h c ij ( τ ) i − . − . . . . c i j ( τ ) P ( c ij ) r i j . . . . h c ij ( τ ) i − . − . . . . c i j ( τ ) P ( c ij ) Figure 5.
Correlations for different time lags, distances and years and their corre-sponding distributions.
The left panels show average of the correlation coefficients (cid:104) c ij ( τ ) (cid:105) asa function of the distance r ij between neighborhoods i and j time lags of τ = 0 weeks, τ = 2weeks and τ = 4 weeks. As depicted, epidemic years (2011, 2012, 2015, and 2016) present a highervalue of (cid:104) c ij ( τ ) (cid:105) when compared to non-epidemic years (2013 and 2014). The right panels show thecorresponding distribution P ( c ij ) of the correlation coefficients c ij ( τ, r ij ). Epidemic years presenta negative skewness for P ( c ij ), while non-epidemic years have no skewness. another more distant region due to human mobility, because mosquitoes at the new locationthen have the possibility of getting in contact with the virus [30–32]. To analyze this effectand identify characteristic correlation length-scales, we study the correlation between the11engue prevalence in neighborhood i and another neighborhood j located within a distanceof r ij . Here, the distance r ij is the distance between the barycenters of neighborhoods i and j extracted from their geographical contours. Specifically, we study the correlations betweenneighborhoods i and j , i.e. c ij ( τ, r ij ) as defined by Eq. (2), for different time lags τ and radii r ij . According to the definition of c ij ( τ, r ij ) in Eq. (2), a correlation length ξ smaller thanour considered system size would lead to an observable decay of c ij ( τ, r ij ) for r ij > ξ . In allconsidered years, the correlations c ij ( τ ) only vary slightly with the distance r ij , as shownin Fig. 5. Within the error bars of our data we do not observe a substantial decay of c ij ( τ )for values of r ij smaller than the system size of about 15 km. We thus conclude that thecharacteristic correlation length ξ is of the order of the system size.We choose three time lags, namely, τ = 0 weeks, τ = 2 weeks and τ = 4 weekswhich are of the order of the transmission time scale of 12 days [1]. The dependence of thecorrelations on different time lags, distances and years is shown in Fig. 5 (left). For no timelag, τ = 0 weeks, or a time lag of τ = 2 weeks, we find that the correlations are substantiallylarger in the epidemic years 2011, 2012, 2015 and 2016 as compared to the non-epidemicyears 2013 and 2014. However, in 2016 the correlations are less pronounced compared toother epidemic years since the outbreaks are widely distributed and their densities rathersmall, as shown in Fig. 3. In the case of the larger time lag, τ = 4 weeks, the averagecorrelation (cid:104) c ij ( τ ) (cid:105) for epidemic years decreases when compared with smaller time lags.This behavior is in accordance with observation, since a time lag of 4 weeks exceeds thedengue transmission period, and we expect to find smaller correlations. In Fig. 5 (right), weshow the corresponding distributions of P ( c ij ) and find that they allow to clearly distinguishbetween epidemic and non-epidemic years for no time lag and a time lag τ = 2 weeks. Inparticular in the case of no time lag, the correlations are substantially larger in epidemicyears compared to the non-epidemic ones.It is unlikely that disease vectors are responsible for such correlation effects over distancesof multiple kilometers due to their limited movement capabilities. In particular, it is knownthat the maximum flight distance reached by the Aedes aegypti from its breeding locationis of the order of approximately 100 meters. On the other hand, humans travel through thedensely populated urban region on a daily basis, and may transfer the virus to Aedes aegyptimosquitoes at different locations. This virus transfer mechanism is therefore compatible withthe large correlation lengths revealed in this study and also agrees well with recent studies12hat suggest that human mobility is a key component of dengue spreading [30–32]. ESTIMATING DISEASE TRANSMISSION PARAMETERS
As mentioned in the previous section, the considered dengue outbreaks in Fortaleza ex-hibit a correlation length ξ which is of the order of the underlying system size. Based onthis observation, we apply mean-field epidemic model, as a first approximation, to furthercharacterize the observed spreading dynamics. In particular, we aim at comparing the dis-ease transmission parameters of dengue outbreaks in Fortaleza with the results of previousstudies [34]. To do so, we consider the two epidemic years 2012 and 2015 and perform aBayesian Markov chain Monte Carlo parameter estimation for a SIR model that has beenpreviously applied in the context of dengue outbreaks in Thailand [34]. More details aboutthe methodology are presented in Appendix A. Such a parameter estimation enables us todetermine the basic reproduction number R of dengue outbreaks in Fortaleza and compareit with the values of other disease outbreaks. The basic reproduction number constitutes animportant epidemiological measure and is defined as the average number of secondary casesoriginating from one infectious individual during the initial outbreak period (i.e., in a fullysusceptible population) [40]. Different models exist to characterize and predict epidemicoutbreaks [40–44].In the case of dengue, some models explicitly incorporate a mosquito population, whereasothers take into account such effects by using an effective spreading rate [34, 45]. Accordingto Ref. [34], an explicit treatment of vector populations just increases the number of modelingparameters and may not lead to a better agreement between model and data. We thereforedo not explicitly model a vector population and consider a SIR model with an effectivespreading rate which has been found to capture essential features of dengue outbreaks inThailand [34]. The governing equations ared S d t = µ H N − β IN S − µ H S, (3)d I d t = β IN S − γ H I − µ H I, (4)d R d t = γ H I − µ H R, (5)where N = S + I + R is the human population size and S = S ( t ), I = I ( t ), R = R ( t )the numbers of susceptible, infected and recovered individuals at time t , respectively. The13uman mortality rate is denoted by µ H and the human recovery rate by γ H . The term µ H N in Eq. (3) corresponds to the birth rate and is chosen such that the population size is keptconstant. Furthermore, the composite human-to-human transmission rate β accounts fortransitions from susceptible to infected. The relation β ≈ mc β H β V µ V (6)connects β with the mosquito-to-human transmission rate per bite β H , human-to-mosquitotransmission rate β V , number of mosquitoes per person m , mean rate of bites per mosquito c , and mosquito mortality rate µ V [34]. The basic reproduction number of the SIR modelis [40] R = βµ H + γ H . (7)For R > R ≤
1. In addition to Eqs. (3)-(5), we define the time evolution of the cumulativenumber of reported dengue cases C = C ( t ) byd C d t = pβ IN S, (8)where p is the fraction of infected individuals that were diagnosed with dengue and reportedto the health officials. We set the population size to N = 2 . µ H = 1 /
76 y − in accordance with the latest World Bank life expectancy estimates [46].A parameter estimation based on Eqs. (3)-(5) and (8) means that we have to determinethe parameters β , γ H , p and the initial fraction of recovered individuals r = R (0) /N thatbest describe the observed dengue outbreaks. More specifically, we use Bayes’ theorem todetermine the posterior parameter distribution P ( θ | D ) ∝ P ( D | θ ) P ( θ ) (9)based on the likelihood function P ( D | θ ) (i.e., the conditional probability of obtaining thedata D for given model parameters θ ) and the prior parameter distribution P ( θ ). Theactual Bayesian Markov chain Monte Carlo algorithm is described in Appendix A.In order to be able to compare our results with the ones of Ref. [34] and to obtain afull infection period, we consider November as our starting month, because the rainfall thentypically reaches a minimum. In Table III, we summarize the inferred model parameters,14edian values, and 95% confidence intervals of the posterior distributions. In addition,we also present the maximum likelihood (ML) estimates which we use to compare the SIRmodel with the actual dengue case data in Fig. 7. We find good agreement between the modelprediction and the actually reported number of dengue cases. Only the dengue outbreakpeak between April and June 2012 is difficult to capture because due to overwhelmingnumber of up to 1000 new infections per day. .
12 0 .
14 0 .
16 0 .
18 0 .
20 0 .
22 0 . β P D F .
10 0 .
12 0 .
14 0 .
16 0 .
18 0 . γ H P D F .
02 0 .
03 0 .
04 0 .
05 0 .
06 0 .
07 0 . p P D F .
00 0 .
05 0 .
10 0 .
15 0 .
20 0 .
25 0 .
30 0 . r P D F Figure 6.
Posterior distributions of SIR model parameters.
The probability density func-tions (PDFs) of the posterior distributions of the modified SIR model as defined by Eqs. (3)-(5)and (8). The two epidemic years 2012 (orange bars) and 2015 (blue bars) have been considered forthe parameter estimation.
After estimating the SIR model parameters for the epidemic years 2012 and 2015, wenow compare the obtained estimates with the ones of Ref. [34] which focuses on dengueoutbreaks in Thailand. In ref [34], the number of reported dengue hemorrhagic fever (DHF)cases, a more severe form of dengue fever, has been used, while we considered the numberof all reported dengue cases. The total number of mentioned DHF cases in Thailand isroughly 70000 with a total population size of 46 . arameter ML Median 95% CI β (cid:0) d − (cid:1) Composite transmission rate 0.1453 (2012)0.1518 (2015) 0.1648 (2012)0.1606 (2015) (0.1409, 0.2138) (2012)(0.1453, 0.1966) (2015) γ H (cid:0) d − (cid:1) Human recovery rate 0.1011 (2012)0.1009 (2015) 0.1121 (2012)0.1064 (2015) (0.1005, 0.1650) (2012)(0.1004, 0.1299) (2015) p Probability of reporting a dengue case 0.0368 (2012)0.0218 (2015) 0.0432 (2012)0.0234 (2015) (0.0348, 0.0609) (2012)(0.0206, 0.0291) (2015) r Initial fraction of humans recovered 0.0594 (2012)0.0618 (2015) 0.0769 (2012)0.0636 (2015) (0.0066, 0.2928) (2012)(0.0045, 0.2609) (2015) R Basic reproduction number 1.4367 (2012)1.5039 (2015) 1.5384 (2012)1.4904 (2015) (1.2184, 1.9584) (2012)(1.3471, 1.7396) (2015)Table III.
Posterior parameter estimations.
Based on given prior SIR parameter distributions,the corresponding posterior distributions for the dengue outbreaks in 2012 and 2015 have been ob-tained using Bayesian Markov chain Monte Carlo sampling. For both years an additional maximumlikelihood (ML) parameter estimate is given. The posterior distributions are characterized by theirmedian values and their 95% confidence intervals (CI). total population size of 2 . p and r are a factor 2-3 larger in our parameter estimation. On the other hand, the MLestimates of β and γ are slightly smaller in Fortaleza. The basic reproduction number asdefined in Eq. (7) (basically the fraction of β and γ H ) is 1 .
44 (2012) and 1 .
50 (2015) andthus exhibits a value within the confidence interval of Ref. [34]. This result suggests thatthe dengue outbreaks in both locations show similar characteristics. In other words, theaverage number of secondary dengue cases originating from one infection is almost identicalin both locations. 16
50 100 150 200 250 300 350 400
Days (November 2011-December 2012) C u m u l a ti v e nu m b e r o f ca s e s r e po r t e d ( t hou s a nd s ) SIRData
Days (November 2014-December 2015) C u m u l a ti v e nu m b e r o f ca s e s r e po r t e d ( t hou s a nd s ) SIRData
Figure 7.
SIR model fit for the epidemic years 2012 and 2015.
The maximum likelihoodfits of the SIR model (black solid lines) as defined by Eqs. (3)-(5) and (8) for the epidemic years2012 (left) and 2015 (right). The cumulative number of reported dengue cases according to ourdata set is indicated by green circles.
DISCUSSION
The re-emergence of dengue and other neglected tropical diseases is a major threat inmany countries. We analyzed the recent dengue outbreaks in Fortaleza as one of the largestBrazilian cities. We identified regions which exhibit a large number of dengue infections andAedes aegypti larvae over different years. Our results show that the characteristic lengthscale of correlation effects between the number of cases at different locations is of the orderof the system size. One possible explanation for this observation is that the human mobility,and not the epidemic vector (the mosquito), is the major mechanism responsible for thedissemination of the virus. Human movement seems to be a crucial factor in the transmissiondynamics, particularly in large tropical cities where successive dengue epidemics have beenrecorded, suggesting that the mosquito is highly adapted, with low intra-urban infestationdifferentials. In addition, we also compared the disease transmission characteristics of dengueoutbreaks in Fortaleza with the ones in Thailand. Our comparison showed that the numberof dengue cases in Fortaleza is much higher compared to the outbreak data of Ref. [34]. Wealso found that the obtained basic reproduction number of dengue outbreaks in Fortalezais within the confidence interval of the value mentioned in Ref. [34]. This means that theaverage number of secondary dengue cases originating from one infection is almost identicalin both locations. Moreover, we classified epidemic and non-epidemic years based on ananalysis of spatio-temporal correlations and their corresponding distributions. Finally, we17ound that the spatial-temporal correlations are of the size of the system likely due to humanmobility.
Appendix A: Bayesian Markov chain Monte Carlo
To compute the parameter distributions of the SIR model defined by Eqs. (3)-(5) and(8), we have to determine the probability distribution P ( θ | D ) of θ = ( β, γ H , p, r ) given ourdata on the cumulative number of reported dengue cases D . According to Bayes’ theorem,we express the posterior distribution P ( θ | D ) ∝ P ( D | θ ) P ( θ ) , (A1)in terms of the likelihood function P ( D | θ ) and the prior parameter distribution P ( θ ). Weassume that the likelihood function P ( D | θ ) is a Gaussian distribution of the error E withzero mean and variance σ , i.e., P ( D | θ ) ∝ exp (cid:18) − E σ (cid:19) . (A2)To compute the error, we discretize the infection period in monthly time intervals describedby the times t i with i ∈ { , , . . . , } and evaluate y i ( θ ) = C ( t i ) /N for a given θ . In thenext step, we compare the prediction of our model y i ( θ ) with the actual cumulative numberof reported dengue cases D i at time t i and use the least square error E = (cid:88) i =1 [ D i − y i ( θ )] (A3)as our error estimate in Eq. (A2) [34]. For a given prior parameter distribution P ( θ ),we compute the posterior distribution P ( D | θ ) using Bayesian Markov chain Monte Carlosampling with a Metropolis update scheme [34, 47]. After initializing the parameter vectorwith θ drawn from the prior distribution P ( θ ), the n th iteration of the algorithm is definedas follows:1. A new parameter set θ ∗ is drawn from the proposal distribution J ( θ ∗ | θ n ).2. The acceptance probability for θ ∗ is computed according to (Metropolis algorithm) r = min (cid:18) P ( θ ∗ | D ) P ( θ n | D ) , (cid:19) = min (cid:18) P ( D | θ ∗ ) P ( θ ∗ ) P ( D | θ n ) P ( θ n ) , (cid:19) . (A4)18. Draw a random number (cid:15) ∼ U (0 ,
1) and set θ n +1 = θ ∗ if (cid:15) < r,θ n otherwise. (A5)This update procedure implies that a new parameter set is always accepted if the newlikelihood function value is greater or equal than the one of the previous iteration, i.e., if P ( D | θ ∗ ) ≥ P ( D | θ n ). For the described Metropolis algorithm, the proposal distribution J ( θ ∗ | θ n ) must be symmetric. According to Ref. [47], a multivariate Gaussian distribution J ( θ ∗ | θ n ) ∼ N (cid:0) θ n | λ Σ (cid:1) (A6)may be used as proposal distribution. The covariance matrix is denoted by Σ and thecorresponding scaling factor by λ . Every 500 iterations, the covariance matrix and thescaling factor are updated using the following update procedure [34, 47]:Σ k +1 = p Σ k + (1 − p )Σ ∗ ,λ k +1 = λ k exp (cid:18) α ∗ − ˆ αk (cid:19) , (A7)where Σ ∗ and α ∗ are the covariance matrix and the acceptance rate of the last 500 iterations,respectively. Furthermore, p = 0 .
25 and the initially Σ = I and λ = 2 . / √ d where I isthe d × d identity matrix and d the number of estimated parameters. The target acceptancerate is [47] ˆ α = .
44 if d = 1 , .
23 otherwise. (A8)For evaluating Eq. (A4), we have to compute the least-square error for the new proposedparameter set θ ∗ according to Eq. (A3) to then obtain the likelihood function value P ( D | θ ∗ )based on Eq. (A2). To obtain good convergence behavior, we set the variance of our like-lihood distribution to σ = 1 /
20. Convergence is measured in terms of the Gelman-Rubintest [47]. Therefore, we consider four independent Markov chains. Every chain is initial-ized with a different random parameter set which is drawn from the corresponding uniformdistributions. The posterior parameter distribution is considered to be converged if the vari-ance between the chains is similar to the variance within the chain (Gelman-Rubin test).19 Iterations − − − − − − E rr o r chain 1chain 2chain 3chain 4 Iterations − − − − − − − E rr o r chain 1chain 2chain 3chain 4 Figure A.1.
Convergence of Bayesian Markov chain Monte Carlo parameter estima-tion.
For both epidemic years 2012 (left) and 2015 (right), four different chains with differentinitial conditions have been considered to perform Bayesian Markov chain Monte Carlo parameterestimation. The algorithm is considered to be converged to the posterior distribution if the variancebetween the chains is similar to the variance within the chains (Gelman-Rubin test) [47].
In Fig. A.1, we illustrate the convergence behavior. After reaching convergence, we pro-duce 10 more samples without updating the covariance matrix and construct the posteriorparameter distribution based on every 5th sample.For both epidemic years, we use four independent Markov chains with different initialparameter sets θ = ( β, γ H , p, r ) that are drawn from uniform posterior distributions. Ifthe variance between the chains is similar to the variance within the chain (Gelman-Rubintest), we consider the posterior distribution to be converged. Every new proposed set ofparameters θ ∗ is accepted with probability (Metropolis algorithm) [47] r = min (cid:18) P ( θ ∗ | D ) P ( θ | D ) , (cid:19) = min (cid:18) P ( D | θ ∗ ) P ( θ ∗ ) P ( D | θ ) P ( θ ) , (cid:19) . (A9)The resulting posterior distributions are shown in Fig. 6. We find similar posterior distri-butions of the parameters β , γ H and r for both epidemic years. The distributions of p (fraction of infected individuals that were diagnosed with dengue and reported to the healthofficials) are different due to the larger number of reported dengue cases in 2012 (38197)compared to 2015 (26176). Declarations vailability of data and material All data generated or analysed during this study are included in this published article[and its supplementary information files].
Competing interests
The authors declare no competing interest.
Acknowledgements
We acknowledge financial support from the Brazilian agencies CNPq, CAPES, the FUN-CAP Cientista Chefe program, and from the National Institute of Science and Technologyfor Complex Systems (INCT-SC) in Brazil. LB acknowledges financial support from theSNF Early Postdoc. Mobility fellowship on “Multispecies interacting stochastic systems inbiology” and the Army Research Office (W911NF-18-1-0345).
Authors’ contributions
All authors designed research. GSS and ASLN collected and prepared the Dengue datasets. SDRS, LB, and JPdCN processed the Dengue data. All authors analyzed and inter-preted the Dengue data. SDSR, LB, ASLN, HH and JSA wrote the manuscript. All authorsread and approved the final manuscript. [1] World Health Organization. Sustaining the drive to overcome the global impact of neglectedtropical diseases; 2013.[2] World Health Organization. Dengue: guidelines for diagnosis, treatment, prevention andcontrol; 2009.[3] Barret ADT, Higgs S. Yellow Fever: A Disease that Has Yet to be Conquered. Annu RevEntomol. 2007;(52):209–229.
4] Charrel RN, et al. Chikungunya outbreaks-the globalization of vectorborne diseases. N EnglJ Med. 2007;(356):769.[5] Brasil P, et al. Zika Virus Infection in Pregnant Women in Rio de Janeiro. N Engl J Med.2016;(375):2321–2334.[6] Hotez PJ, Bottazzi ME, Franco-Paredes C, Ault SK, Periago MR. The Neglected TropicalDiseases of Latin America and the Caribbean: A Review of Disease Burden and Distributionand a Roadmap for Control and Elimination. PLoS Negl Trop Dis. 2008;(2):e300.[7] World Health Organization. A global brief on vector-borne diseases; 2014.[8] Enslen AW, Lima Neto AS, Castro MC. Infestation measured by Aedes aegypti larval surveysas an indication of future dengue epidemics: an evaluation for Brazil. Transactions of TheRoyal Society of Tropical Medicine and Hygiene. 2020;.[9] MacCormack-Gelles B, Neto ASL, Sousa GS, do Nascimento OJ, Castro MC. Evaluation ofthe usefulness of Aedes aegypti rapid larval surveys to anticipate seasonal dengue transmissionbetween 2012–2015 in Fortaleza, Brazil. Acta Tropica. 2020;205:105391.[10] Sanchez L, Vanlerberghe V, Alfonso L, del Carmen Marquetti M, Guzman MG, Bisset J,et al. Aedes aegypti larval indices and risk for dengue epidemics. Emerging infectious diseases.2006;12(5):800.[11] Cromwell EA, Stoddard ST, Barker CM, Van Rie A, Messer WB, Meshnick SR, et al. Therelationship between entomological indicators of Aedes aegypti abundance and dengue virusinfection. PLoS neglected tropical diseases. 2017;11(3).[12] Focks DA, Chadee DD. Pupal survey: an epidemiologically significant surveillance methodfor Aedes aegypti: an example using data from Trinidad. The American journal of tropicalmedicine and hygiene. 1997;56(2):159–167.[13] Focks DA, et al. A review of entomological sampling methods and indicators for denguevectors. World Health Organization; 2004.[14] Gubler DJ. Dengue and dengue hemorrhagic fever. Clinical microbiology reviews.1998;11(3):480–496.[15] Bowman LR, Runge-Ranzinger S, McCall P. Assessing the relationship between vector in-dices and dengue transmission: a systematic review of the evidence. PLoS neglected tropicaldiseases. 2014;8(5).[16] MacCormack-Gelles B, Neto ASL, Sousa GS, Nascimento OJ, Machado MM, Wilson ME, t al. Epidemiological characteristics and determinants of dengue transmission during epidemicand non-epidemic years in Fortaleza, Brazil: 2011-2015. PLoS neglected tropical diseases.2018;12(12):e0006990.[17] Tun-Lin W, Kay B, Barnes A, Forsyth S. Critical examination of Aedes aegypti indices:correlations with abundance. The American journal of tropical medicine and hygiene.1996;54(5):543–547.[18] Bouzid M, Brainard J, Hooper L, Hunter PR. Public health interventions for Aedes controlin the time of Zikavirus–a meta-review on effectiveness of vector control strategies. PLoSneglected tropical diseases. 2016;10(12).[19] Stoddard ST, Forshey BM, Morrison AC, Paz-Soldan VA, Vazquez-Prokopec GM, Astete H,et al. House-to-house human movement drives dengue virus transmission. Proceedings of theNational Academy of Sciences. 2013;110(3):994–999.[20] Vazquez-Prokopec GM, Bisanzio D, Stoddard ST, Paz-Soldan V, Morrison AC, Elder JP, et al.Using GPS technology to quantify human mobility, dynamic contacts and infectious diseasedynamics in a resource-poor urban environment. PloS one. 2013;8(4).[21] Vazquez-Prokopec GM, Kitron U, Montgomery B, Horne P, Ritchie SA. Quantifying thespatial dimension of dengue virus epidemic spread within a tropical urban environment. PLoSneglected tropical diseases. 2010;4(12).[22] Vazquez-Prokopec GM, Stoddard ST, Paz-Soldan V, Morrison AC, Elder JP, Kochel TJ,et al. Usefulness of commercially available GPS data-loggers for tracking human movementand exposure to dengue virus. International journal of health geographics. 2009;8(1):68.[23] Organization PAH. Technical document for the implementation of interventions based ongeneric operational scenarios for Aedes aegypti control. PAHO Washington, DC; 2019.[24] Organization WH, et al. Handbook for integrated vector management. World Health Orga-nization; 2012.[25] Quintero J, Brochero H, Manrique-Saide P, Barrera-P´erez M, Basso C, Romero S, et al.Ecological, biological and social dimensions of dengue vector breeding in five urban settingsof Latin America: a multi-country study. BMC infectious diseases. 2014;14(1):38.[26] Caprara A, Lima JWdO, Marinho ACP, Calvasina PG, Landim LP, Sommerfeld J. Irregularwater supply, household usage and dengue: a bio-social study in the Brazilian Northeast.Cadernos de saude publica. 2009;25:S125–S136.
27] Guzzetta G, Marques-Toledo CA, Ros`a R, Teixeira M, Merler S. Quantifying the spatial spreadof dengue in a non-endemic Brazilian metropolis via transmission chain reconstruction. NatCommun. 2018;9(1):2837.[28] Antonio FJ, Itami AS, de Picoli S, Teixeira JJV, dos Santos Mendes R. Spatial patterns ofdengue cases in Brazil. PloS one. 2017;12(7):e0180715.[29] Stolerman L, Maia P, Kutz JN. Data-Driven Forecast of Dengue Outbreaks in Brazil: ACritical Assessment of Climate Conditions for Different Capitals. arXiv:170100166. 2016;.[30] Stoddard ST, et al. House-to-house human movement drives dengue virus transmission. ProcNatl Acad Sci. 2013;(110):994–999.[31] Jr RCR, Stoddard ST, Scott TW. Socially structured human movement shapes dengue trans-mission despite the diffusive effect of mosquito dispersal. Epidemics. 2014;(6):30–36.[32] Chris M Stone DMF Samantha R Schwab, Fefferman NH. Human movement, cooper-ation and the effectiveness of coordinated vector control strategies. J R Soc Interface.2017;(14):20170336.[33] Lindoso JAL, Lindoso AABP. Neglected tropical diseases in Brazil. Rev Inst Med trop SPaulo. 2009;(51):247–253.[34] Pandey A, Mubayi A, Medlock J. Comparing vector–host and SIR models for dengue trans-mission. Math Biosci. 2013;246(2):252–259.[35] Fortaleza Daily Disease Monitoring System (SIMDA). https://pt.climate-data.org/america-do-sul/brasil/ceara/fortaleza-2031/, retrieved June 3, 2020;.[36] Minist´erio da Sa´ude do Brasil. Diretrizes nacionais para preven¸c˜ao e controle de epidemias dedengue. 2009;.[37] https://cidades.ibge.gov.br/brasil/ce/fortaleza/panorama, retrieved December 9, 2018;.[38] https://python-visualization.github.io/folium/, retrieved May 4, 2020;.[39] https://pt.climate-data.org/america-do-sul/brasil/ceara/fortaleza-2031/, retrieved December9, 2018;.[40] Keeling MJ, Rohani P. Modeling Infectious Diseases in Humans and Animals. PrincetonUniversity Press; 2008.[41] B¨ottcher L, Woolley-Meza O, Ara´ujo NAM, Herrmann HJ, Helbing D. Disease-induced re-source constraints can trigger explosive epidemics. Sci Rep. 2015;5:16571.[42] B¨ottcher L, Woolley-Meza O, Goles E, Helbing D, Herrmann HJ. Connectivity disruption parks explosive epidemic spreading. Phys Rev E. 2016;93(4):042315.[43] B¨ottcher L, Nagler J, Herrmann HJ. Critical Behaviors in Contagion Dynamics. Phys RevLett. 2017;118:088301.[44] B¨ottcher L, Andrade J, Herrmann HJ. Targeted Recovery as an Effective Strategy againstEpidemic Spreading. Sci Rep. 2017;7(1):14356.[45] Fitzgibbon W, Morgan J, Webb G. An outbreak vector-host epidemic model with spatial struc-ture: the 2015–2016 Zika outbreak in Rio De Janeiro. Theor Biol Med Model. 2017;14(1):7.[46] https://data.worldbank.org/indicator/SP.DYN.LE00.IN, retrieved December 9, 2018;.[47] Gelman A, Stern HS, Carlin JB, Dunson DB, Vehtari A, Rubin DB. Bayesian data analysis.Chapman and Hall/CRC; 2013.parks explosive epidemic spreading. Phys Rev E. 2016;93(4):042315.[43] B¨ottcher L, Nagler J, Herrmann HJ. Critical Behaviors in Contagion Dynamics. Phys RevLett. 2017;118:088301.[44] B¨ottcher L, Andrade J, Herrmann HJ. Targeted Recovery as an Effective Strategy againstEpidemic Spreading. Sci Rep. 2017;7(1):14356.[45] Fitzgibbon W, Morgan J, Webb G. An outbreak vector-host epidemic model with spatial struc-ture: the 2015–2016 Zika outbreak in Rio De Janeiro. Theor Biol Med Model. 2017;14(1):7.[46] https://data.worldbank.org/indicator/SP.DYN.LE00.IN, retrieved December 9, 2018;.[47] Gelman A, Stern HS, Carlin JB, Dunson DB, Vehtari A, Rubin DB. Bayesian data analysis.Chapman and Hall/CRC; 2013.