Impact of human mobility on the emergence of dengue epidemics in Pakistan

Edited by Burton H. Singer, University of Florida, Gainesville, FL, and approved August 6, 2015 (received for review April 2, 2015)
September 8, 2015
112 (38) 11887-11892

Significance

Dengue virus has rapidly spread into new human populations due to human travel and changing suitability for the mosquito vector, causing severe febrile illness and significant mortality. Accurate predictive models identifying changing vulnerability to dengue outbreaks are necessary for epidemic preparedness and containment of the virus. Here we show that an epidemiological model of dengue transmission in travelers, based on mobility data from ∼40 million mobile phone subscribers and climatic information, predicts the geographic spread and timing of epidemics throughout the country. We generate fine-scale dynamic risk maps with direct application to dengue containment and epidemic preparedness.

Abstract

The recent emergence of dengue viruses into new susceptible human populations throughout Asia and the Middle East, driven in part by human travel on both local and global scales, represents a significant global health risk, particularly in areas with changing climatic suitability for the mosquito vector. In Pakistan, dengue has been endemic for decades in the southern port city of Karachi, but large epidemics in the northeast have emerged only since 2011. Pakistan is therefore representative of many countries on the verge of countrywide endemic dengue transmission, where prevention, surveillance, and preparedness are key priorities in previously dengue-free regions. We analyze spatially explicit dengue case data from a large outbreak in Pakistan in 2013 and compare the dynamics of the epidemic to an epidemiological model of dengue virus transmission based on climate and mobility data from ∼40 million mobile phone subscribers. We find that mobile phone-based mobility estimates predict the geographic spread and timing of epidemics in both recently epidemic and emerging locations. We combine transmission suitability maps with estimates of seasonal dengue virus importation to generate fine-scale dynamic risk maps with direct application to dengue containment and epidemic preparedness.
Dengue is the most rapidly spreading mosquito-borne disease worldwide (1, 2). Half the global population now lives in at-risk regions for dengue virus transmission, due to the wide distribution of the mosquito vector, Aedes aegypti, which thrives in peri-urban areas and transmits the virus between humans (3). Dengue virus can cause acute febrile illness and carries the risk of severe disease, hospitalization, and shock syndrome, especially in clinical settings with little experience treating dengue patients. There is currently no specific therapeutic protocol for, or vaccine against, infection (1). Current control measures focus on vector control, although these measures are often logistically difficult and have shown varying efficacy in controlling epidemics (4). In the absence of effective prevention and treatment, public health system preparedness remains the single most important tool for minimizing morbidity and mortality as dengue epidemics spread beyond endemic areas (5, 6).
The introduction of dengue into new populations is mediated by travel of infected individuals to areas that can support transmission, because mosquito vectors move only short distances during their lifespans (3, 712). International travel to endemic countries has resulted in imported cases and outbreaks in Europe and the Americas (2, 8, 10, 13). Local variation in transmission, within a single city for example, is also driven by mobility patterns of individuals on short timescales (7). Forecasting methods are needed to spatially target interventions and epidemic preparedness measures that reflect both the changing temporal risks of importation and environmental suitability that go beyond solely climate-based methods (14).
Dengue has long been endemic in most Southeast Asian countries (1), but has more recently emerged in parts of the Middle East and South Asia, including Pakistan (15, 16). In Pakistan, the transmission of dengue viruses was largely confined to the southern city of Karachi until 2011 when a large dengue epidemic with over 20,000 cases occurred in the northeastern city of Lahore (16), causing significant morbidity and mortality. In 2013, a second large epidemic occurred in northeastern Pakistan in Punjab and Khyber-Pakhtunkhwa (KP) provinces, establishing the region as an emerging focus of seasonal dengue epidemics. It has been hypothesized that the recent geographic expansion of A. aegypti mosquito vectors, changing environmental suitability, and human importation of dengue from endemic regions all contributed to the emergence of dengue in northern areas (17). Pakistan is therefore representative of many countries that are on the verge of countrywide endemic dengue transmission and are struggling to contain its emergence into previously dengue-free regions.
Measuring changing risks of importation events that spark epidemics has been extremely challenging on the refined temporal and spatial scales necessary to inform local policies (18). Being able to predict when to prepare surveillance systems and health facilities for dengue outbreaks could dramatically reduce the morbidity and mortality associated with epidemics and would allow policy makers to pinpoint regions that are particularly vulnerable to imported cases, for vector control. Mobile phone data offer direct measures of human aggregation and movement and represent a unique source of information on the human determinants of the geographic expansion of emerging epidemic diseases like dengue. Here, we conduct a retrospective epidemiological analysis of large dengue outbreaks in Pakistan in 2013, to examine the predictive ability of an epidemiological model that integrates human mobility from the largest mobile phone dataset analyzed to date with climate information. We show that within-country human mobility predicts emerging epidemics in Pakistan, and epidemiological models incorporating this type of data can predict the spatial extent and timing of outbreaks, providing a new approach to forecasting.

Results

Human Mobility in Pakistan Does Not Conform to Standard Model Predictions.

To measure human travel patterns underlying the spread of dengue virus across Pakistan in 2013, we estimated the mobility of 39,785,786 mobile phone subscribers between geolocated mobile phone towers in Pakistan between June 1 and December 31 of 2013 (representing ∼22% of the population; Materials and Methods). Daily locations and movements were aggregated to measure travel between 356 small, politically defined areas called tehsils (Fig. 1A and Materials and Methods). We compared our data to gravity models of mobility, developed from transportation theory and commonly used to parameterize infectious disease frameworks, to assess whether observed mobility measured using mobile phone data significantly improves upon this standard approach. We have focused on intracountry mobility patterns (Discussion).
Fig. 1.
Human mobility dynamics in Pakistan. (A) Population density (red, high density; yellow, low density) and mobile phone tower coverage from the mobile phone operator in Pakistan (colored in gray) per tehsil. (B) The top routes of travel between pairs of tehsils in Pakistan. A line is drawn if at least 20,000 trips occurred between the origin and destination between June and December 2013. The top routes occur between Karachi and cities in northern Punjab province, particularly Lahore tehsil. (C) Relative direction and volume of travel. For each trip, we calculated the distance traveled from the origin and the destination. The origin location was centered at 0,0 and the longitude distance and latitude distance to the destination are shown. Although many trips occurred over short distances, a substantial amount of travel occurred between the southeastern and northern parts of the country, reflecting the geography and population distribution of Pakistan.
The sampled population was extremely mobile. We estimated that between 2.4 million and 4.8 million subscribers traveled between tehsils each day (95% quantile interval: 3.1–4.6 million; details in Materials and Methods). Most travel followed a NW–SE corridor along the major highways (Fig. 1B). Large volumes of travel occurred to and from Karachi, a major population and economic hub of Pakistan, with ∼710,000 subscribers traveling to or from the city each day on average (95% quantile interval: 570,000–813,000; Fig. S1). In contrast to expectations of standard mobility models, there was almost no decay in travel with increasing distance (correlation coefficient: −0.064, P < 0.001), although the most frequent destinations for travel were often in a nearby tehsil (Fig. 1C). This pattern reflects the topography, road infrastructure, and population distribution in Pakistan, with the largest cities outside Karachi being located in the northern part of Punjab province that includes the Rawalpindi/Islamabad metropolitan area, Lahore, and Faisalabad. Although there was a decrease in overall movement during Ramadan, we did not observe a systematic difference in the amount or direction of travel patterns in the weeks before and after the holiday (Fig. S1).
Fig. S1.
The mobility values derived from the mobile phone data. (A) The mobility estimates for travel from Karachi to Lahore. Using the mobile phone data, we calculated both the number of trips between tehsils (A, Top Left) and the daily number of active subscribers: flux value (A, Top Right). The number of trips was normalized by the daily number of active subscribers (A, Bottom Left) that was used to estimate the importation of infected travelers from Karachi. Our dataset runs from June to December 2013 and we inferred the mobility estimates (gray) for January–June of that year (Materials and Methods) (A, Bottom Right). (B) The most frequently traveled routes based on the mobile phone data. We compared the most traveled routes (top 10,000 routes) 1 mo before Ramadan (B, Left) (June), during Ramadan (July) (B, Center), and after Ramadan (B, Right) (September). The top routes remain consistent throughout the major holiday, in particular the large amount of travel between Karachi and Punjab province. A line is drawn if more than 10,000 trips were taken over a route. (C) Country-wide population flux values per day. Total population flux per day was measured (red line). Ramadan is shown in gray. We have adjusted the mobile phone data (CDRs) (black line shows adjustment) to account for a sharp drop in the raw CDRs after day 145 and fitted the adjusted data [spline, blue line; moving average (MA) trend, purple line], using two common time series methods (a MA and smoothing spline) to smooth the data. For the remaining analysis, we use these trend lines from these time series smoothing methods.

Epidemiological Modeling of Dengue Epidemics in 2013.

There were 15,535 reported dengue cases in 82 tehsils over 7 mo in Pakistan in 2013 (Fig. 2A and Materials and Methods). Peak timing of the epidemic varied by location, and the majority of cases occurred in and around Karachi, in the northern district of Swat, and in the cities of Lahore and Rawalpindi (Fig. 2A and Table S1). About half of the dengue cases reported occurred in the Mingora area of Swat (KP province), marking the first major outbreak in the region [n = 7,950, compared with a previous maximum of 300 cases reported in KP in 2011 (16)].
Fig. 2.
Dengue epidemiology in Pakistan in 2013. (A) The location of dengue cases throughout Pakistan and the number of cases per week by tehsil. Tehsils that reported at least 15 cases are shown on the map with corresponding color shown in the time series. The majority of cases were reported in Karachi (gray), Lahore (blue), and Mingora (orange). The dengue season in the entire country lasted 35 wk, with the first reported case in Karachi during week 18 (end of April). (B) The reported cases (red), temperature (blue), and model fit (black) for Karachi are shown. Using the case and temperature data, the human and vector population dynamics were modeled (Materials and Methods).
Table S1.
The timing and number of reported dengue cases per tehsil in 2013 used in the analysis
TehsilDate of first case in 2013Date of peak cases in 2013Date of last case in 2013Total no. cases in 2013
BabuzaiAugust 19September 17November 184,029
LahoreMay 13November 16December 261,538
RawalpindiAugust 30November 4December 19883
Matta SebujniAugust 22September 15November 12864
BarikotAugust 20September 2November 15731
Matta KhararaiAugust 20September 18November 12678
KhwazakhelaAugust 21September 6November 17413
ChiksarAugust 22September 21October 5386
KabalAugust 23September 8November 5357
AlpuriAugust 22September 18October 8161
KalamAugust 24September 8September 3060
DassuSeptember 3September 17November 956
CharbaghAugust 22September 10October 1049
FerozewalaAugust 5August 5December 1346
AlaiOctober 5October 6November 137
BahrainSeptember 3September 21November 434
SheikhupuraAugust 7September 12November 2724
MamundSeptember 1September 13October 2522
FaisalabadApril 30November 9December 2020
PeshawarOctober 5October 5October 2218
MuridkeAugust 21September 18December 114
BahawalpurSeptember 19November 9November 2713
KasurSeptember 25November 13November 3011
KatlangSeptember 12September 12October 1411
CharsaddaSeptember 1September 5October 259
MultanSeptember 5October 8December 199
BalakotAugust 27September 8September 278
BishamAugust 24September 5September 218
MartungSeptember 6September 17October 248
GujranwalaSeptember 5November 13December 56
ChagharzaiAugust 30September 17September 175
SahiwalSeptember 11November 1November 215
ChakwalOctober 12October 12November 74
NarowalNovember 4November 4November 294
SialkotAugust 20December 12December 124
T.T.SinghJune 29September 23November 74
AttockOctober 10October 24October 243
BahawalnagarSeptember 20September 20November 93
D.G.KhanSeptember 19December 14December 143
DepalpurSeptember 15September 15November 233
HasilpurOctober 4October 16October 163
KahutaOctober 5October 5November 213
Rahimyar KhanSeptember 16September 16October 313
SargodhaSeptember 20December 3December 33
Arif WalaSeptember 2October 27October 272
BurewalaOctober 21December 14December 142
ChiniotSeptember 17September 17October 222
ChunianOctober 29October 29November 162
DirOctober 3October 5October 52
Karor Lal EsanAugust 31September 16September 162
LodhranOctober 14October 14November 192
PakpattanAugust 31August 31September 242
ShakargarhAugust 30August 30December 232
TandlianwalaAugust 5August 5October 92
VehariOctober 31November 1November 12
BhalwalNovember 1November 1November 11
ChichawatniNovember 28November 28November 281
GujratAugust 21August 21August 211
HafizabadOctober 4October 4October 41
HaroonabadNovember 17November 17November 171
JahanianDecember 5December 5December 51
JatoiSeptember 19September 19September 191
KhanewalNovember 3November 3November 31
KharianNovember 2November 2November 21
KhushabNovember 13November 13November 131
Kot AdduAugust 26August 26August 261
LayyahNovember 14November 14November 141
MailsiJuly 16July 16July 161
MardanNovember 1November 1November 11
Mian ChannuOctober 31October 31October 311
Nankana SahibAugust 5August 5August 51
PasrurOctober 29October 29October 291
Pindi BhattianAugust 6August 6August 61
SafdarabadOctober 4October 4October 41
SillanwaliNovember 19November 19November 191
TimergaraOctober 31October 31October 311
YaOctober 22October 22October 221
We first fitted an ento-epidemiological model to the reported dengue cases in Southern Pakistan, where transmission occurs year round, and to case data in Lahore and Swat in northern Pakistan, where transmission is seasonal due to climatic variation. Southern Pakistan has year-round climatic suitability for dengue vectors and is therefore the most likely source of exported cases to other parts of the country, where greater seasonal temperature extremes limit suitability (16) (SI Text). Although importation from international regions is technically a possibility, given the fairly restrictive political borders we assume here that importation will be negligible compared with the within-country importation rate. We fitted an ento-epidemiological model of dengue dynamics in Karachi (Fig. 2B), with vector dynamics being determined by temperature (19). We fitted the mosquito-biting rate (a=0.66) and reporting rate (SI Text), yielding parameter values that are consistent with endemic dengue virus transmission (Materials and Methods and SI Text).
In northern Pakistan, we used our epidemiological framework to estimate the timing of the introduction of the first case that sparked the epidemics in each area, to compare these estimates with our mechanistic model of imported infections using mobile phone data. We performed a sensitivity analysis to estimate an expected range of dates for the introduction of dengue to Lahore and Mingora, where most cases were concentrated (Materials and Methods and Fig. 3 A and B). Interestingly, Lahore had suffered its first major outbreak 2 y before this epidemic, and we expect that immunity may have played a significant role in determining transmission dynamics (16). Mingora, on the other hand, represents an effectively naive population. In the absence of serotype information, we took the simplest approach and assumed each population was immunologically naive; however, we hypothesize that the delay between the first cases and the peak of the epidemic in Lahore may have been caused by immunity from the 2011 outbreak. We estimated that the first case was introduced to Lahore during the second week of May (between days 124 and 130; SI Text and Fig. 3A) a few days earlier than the first reported case (day 133). In Mingora, on the other hand, we estimated that the first introductions likely occurred in August (between days 202 and 231) (Fig. 3B), a few weeks before the first reported case in the city.
Fig. 3.
Mobility estimates derived from mobile phone data predict the timing of introduced cases around the country that spark epidemics. (A and B) The estimated introduced cases from Karachi to (A) Lahore (total dengue cases: 1,538) and (B) Mingora (total dengue cases: 4,029). The estimated introductions (assuming 30% of individuals travel, a 2% reporting rate, and a probability of 0.01) from the mobile phone data (boxplot in blue), the diffusion model (boxplot in green), actual case data (red), and estimated dengue suitability (gray) are shown. Dengue suitability was defined based on temperature and relative humidity, using a measure that is linearly proportional to vectoral capacity (Materials and Methods). Values near zero are unsuitable for dengue transmission. For Lahore and Mingora, the estimated introduction from the case data alone is shown (red cross-hatched box). In all instances, the mobile phone data were able to predict the timing of the first introduced case in each tehsil. An arrow indicates the week of the first reported case in each tehsil.

Models of Dengue Virus Importation Based on Mobile Phone Data Accurately Predict the Spatial Extent and Timing of Epidemics.

We next modeled the number of infected individuals traveling from the endemic areas in southern Pakistan to all other tehsils, using different approaches to characterize mobility: direct observations from the mobile phone data and various modified gravity models of travel (Materials and Methods) (10, 19). To compare the performance of the mobile phone data against the next best alternative, we used a parameter-free gravity model (referred to as the diffusion model) that is equivalent to a population-weighted spatial diffusion model (Materials and Methods) based on the travel time distance between the origin and the destination. In addition to these models, we fitted a gravity model to the mobile phone data, to determine whether simple adjustments would significantly improve our predictions (SI Text and Fig. S2).
Fig. S2.
The various distance measures and gravity model fits to the mobile phone data. We calculated the Euclidean distance (kilometers), road distance (gray, in kilometers), and travel time distance (red, in minutes) between tehsil centroids. (A) The three measures were highly correlated (Pearson’s correlation coefficient between Euclidean distance and road distance, 0.988; between Euclidean distance and travel time distance, 0.99). We then fitted gravity models to the mobile phone data, using each distance measure: (B) Euclidean distance, (C) road distance, and (D) travel time distance. Shown is the normalized predicted amount of travel from each gravity model compared to the mobile phone data.
We compared the timing of predicted importation events from endemic areas in southern Pakistan, based on our mobility model, to the estimated first dengue case inferred from case report data from Lahore and Mingora (SI Text). The timing of the first dengue case estimated from the epidemiological data overlapped well with the predicted introductions from southern Pakistan to Lahore from the importation model (Fig. 3A), with the first introductions occurring approximately 1 mo earlier. In Mingora, the predicted timing of imported infections using our model occurred 2 wk before the first reported case, consistent with the serial interval for dengue and the estimates of the first case from epidemiological data (Fig. 3B). Crucially, the diffusion model does not predict any introductions from endemic areas in southern Pakistan to Mingora. Thus, travel patterns measured using mobile phone data predict introduction events consistent with outbreaks in both emerging (Lahore) and previously dengue-free (Mingora) regions. Our ability to measure these importations in more remote places like Mingora was somewhat sensitive to the Karachi model fit, in particular the reporting rate, although the mobile phone data are always able to predict earlier, more frequent, and more accurate introductions than the diffusion model (SI Text). The modeled interaction between seasonal variations in vectorial capacity and the dynamics of importation events provide accurate predictions about the location and timing of epidemics in different epidemiological settings and regions of the country.
We combined our estimates of imported cases with an index of climatic suitability for dengue vectors (19) across Pakistan to predict the potential for spread of the virus in areas with introduced cases and assessed how well the spatiotemporal dynamics of the epidemic were captured. In general, the spatial extent and epidemic timing of cases were predicted more accurately using mobile phone data (Fig. 4) compared with the diffusion models, which both predicted early spread along the major highway connecting Karachi to other highly populated areas.
Fig. 4.
The estimated spatial spread of imported dengue. Using the modeled dengue dynamics in Karachi and mobility measured from the mobile phone data or a diffusion model, we estimated the time of the first introduced case to the rest of the country. The mobile phone data predict the earliest introductions in eastern Pakistan near Lahore and inland toward Swat Valley (Mingora). In comparison, the mobility model predicts early introductions in southern Pakistan with few introductions in Mingora. These differences are highlighted in the difference in predictions plot—the number of days earlier (red) from the mobile phone predictions or earlier (yellow) from the diffusion model (without the mobile phone data).

A New Approach to Dynamic Risk Mapping for Epidemic Preparedness.

We next constructed risk maps to identify areas of the country that were vulnerable to epidemics due to a combination of seasonally varying climatic suitability based on temperature and relative humidity (vectorial capacity index) for dengue transmission and seasonal patterns of importation of dengue virus from southern Pakistan. By combining the daily suitability index (Fig. 5A) with estimated introductions of infected travelers from endemic areas we created a composite epidemic risk measure that varied in space and time and showed a distinctively different pattern of risk than a climate-based suitability map alone (Fig. 5B). In southern Pakistan outside of Karachi, for example, a low risk of epidemic dengue results from low numbers of infected travelers visiting the area. Along the Afghanistan border, on the other hand, low risk of epidemics is due to temperature extremes that prevent continuous dengue transmission. The highest-risk regions are in Punjab and KP provinces in the northeast of the country, where a combination of frequent introductions from Karachi and long periods of sustained climatic suitability for transmission create the necessary conditions for outbreaks.
Fig. 5.
Dynamic risk mapping for dengue epidemics in Pakistan. (A) Average dengue suitability in Pakistan during 2013 based on temperature and humidity. The southern part of the country is the most suitable for dengue (high values are shown in red). (B) Epidemic risk, measured as a composite of vector suitability and frequency of introductions from endemic areas in southern Pakistan. Places in red are the most at risk whereas those in blue have low risk of epidemics. (C) Variation in the timing of dengue risk by location (per tehsil). A box is drawn if suitability multiplied by estimated introduction events is greater than 0, indicating a nonnegative risk. Each line shown is a tehsil, with the y axis corresponding to the ranked distance from Karachi. The dengue risk is shown using the mobile phone data (blue) or the diffusion model (green) to estimate introduction events from Karachi. In general, the mobile phone data predict earlier risk to tehsils farther away from Karachi, in particular tehsils in Punjab (Faisalabad and Lahore, for example), than the diffusion model data. The diffusion model predicts the earliest risk to nearby tehsils such as Malir.
The timing of epidemic risk varied considerably in different tehsils (Fig. 5C), with the risk of outbreaks varying by up to 4 mo in different regions. Based on the mobile phone data, epidemic risk occurred much earlier in tehsils in northern regions relatively far from Karachi, in contrast to the diffusion model. Because preparation of local health facilities and surveillance for dengue cases are key components of national response to outbreaks, this variation in timing has important implications for decision making. In particular, given the 2-wk serial interval of dengue (20), mobile phone and disease surveillance data could be used to produce dynamic risk maps and provide location-specific predictions about future risk of epidemics in near real time, giving policy makers time to prepare.

Discussion

For emerging epidemic diseases like dengue, changing environments and dynamic population movements create conditions for geographic expansion of pathogens into naive populations. Predicting the changing patterns of risk due to these emerging public health threats will be critical to containing their spread and preparing for outbreaks. Mobile phone data provide dynamic population mobility estimates that can be combined with infectious disease surveillance data and seasonally varying environmental data to map these changing patterns of vulnerability in a country where dengue outbreaks are emerging and irregular in many regions. Because these data are continuously being collected by mobile phone operators, these methods could be integrated into national control programs in near real time.
A major limitation of mobile phone data generated by national operators is the difficulty in capturing cross-border travel patterns. International travel may be less relevant for dengue dynamics in Pakistan than in other dengue-endemic countries, however. Less than 2% of the entire Pakistan population are migrants (where their country of origin is not Pakistan) (21). More than half of these migrants are from Afghanistan (2.3 million) (21). Although Pakistan borders four countries (India, Afghanistan, China, and Iran), all aside from India are not endemic for dengue near the Pakistan borders, making importations of dengue cases to Pakistan unlikely (3, 2224). Estimates of travel between India and Pakistan are difficult to obtain, but travel restrictions between the countries and a limited number of international flights may make travel difficult. Another potential issue with this type of data is the sampling frame, because mobile phone ownership may not represent a representative sample of the population, although we have previously found in an African context that mobility estimates from this kind of data were not significantly affected by ownership bias (25). We expect that the substantial market share and widespread coverage of the mobile phone operator make the possibility of systematic skew in our estimates unlikely, but this cannot be ruled out.
Uncertainties are inherent to any epidemiological model of dengue transmission, in particular regarding the large asymptomatic reservoir and the multiple circulating serotypes of the virus, which impact the ability to estimate both the fraction of the population that is susceptible and the reporting rate accurately (26, 27). In the absence of serotype information, for example, we chose not to include immunological covariates in our model. As Pakistan moves toward an endemic transmission setting in the northeast of the country, however, serotype data will be increasingly important to determine how strain dynamics are likely to affect the likelihood of outbreaks (28). Indeed, we believe that the substantial lag time in Lahore we observed was probably due to the recent large outbreak in 2011, which may have generated substantial immunity and delayed the 2013 epidemic (16, 29). These uncertainties will continue to make forecasting challenging, although serotype data could improve model accuracy.
Given these difficulties, forecasting dengue epidemics will always need to encompass substantial stochastic variation, despite its importance for targeting limited public health resources. Currently, vector control programs in Pakistan begin during the monsoon season uniformly across the country. We believe that the large estimated lead times this approach offers could aid control managers, providing an early warning system. This approach provides policy-relevant, real-time information about where and when to expect dengue epidemics and therefore how to effectively target interventions, surveillance, and clinical response.

Materials and Methods

Population Data.

Pakistan’s large population (182 million) is broadly divided into one capital territory, Federally Administered Tribal Area (FATA), Gilgit-Baltistan, Azad Jammu and Kashmir, and four provinces, which are further subdivided into 388 tehsils (equivalent to administrative unit 3 from 2008; Fig. 1A). Karachi is the most populated city in the country and located along the southern coast. The majority of the population in southern Pakistan lives along the Indus river, whereas in the northern half the population lives in an arc between Faisalabad and Peshawar that includes the major population centers of Lahore, Islamabad, and Rawalpindi. We used population data obtained from worldpop.org.uk.

Dengue Data.

Data for Punjab province were collected by the Provincial Health Department, whereas the District Health Office collected case data for Swat District (Table S1) daily. All public and private hospitals, health clinics, and laboratories reported any case of a patient presenting with dengue symptoms. The data listed each patient who presented to a hospital or clinic, regardless of whether the patient was admitted in the province or district, seeking treatment for high fever, body aches, petechiae, and low platelet counts (with a cutoff value for thrombocytopenia of <50,000/mm3). These cases were then confirmed by using IgM or NS 1-Ag enzyme-linked immunosorbent assay (ELISA). The data reported from Karachi were based on official public releases from the Sindh Health department. All case data were deidentified and aggregated to the tehsil level.

Mobile Phone Data.

We analyzed all voice-based, originated, call data records (CDRs) from 39,785,786 subscriber SIMs (security information management) over a 7-mo period, from June 1 to December 31, 2013. The mobile operator has the largest coverage of tehsil headquarters (352 in the dataset of 388 total tehsils) across Pakistan, particularly in rural areas (two-thirds of the network). To comply with national laws and regulations of Pakistan and the privacy policy of the Telenor group, the following measures were implemented to preserve the privacy rights of Telenor Pakistan’s customers: (i) The CDR/mobility data were processed on a backup and recovery server made available by Telenor Pakistan. Only Telenor employees have access to the detailed CDR/mobility data. (ii) Given the server arrangements, no detailed CDR/mobility data were taken out of Pakistan or left the premises of Telenor Pakistan. (iii) The processing of the detailed CDR/mobility data resulted in aggregations of the data on a tower-level granularity that was accessed only by Telenor employees. Further spatial aggregations to the level of the tehsil were made available to the remaining coauthors.
On average, 28 million subscriber SIMs were recorded as active on a given day, and of these 15.2 million subscribers generated outgoing, voice-event CDRs that encoded location information. At the time of data acquisition, the mobile phone operator had approximately a 25% market share (22% of the population) and was the second largest provider of mobile telecommunication services in Pakistan. Multi-SIM activity is common in Pakistan, but we believe that this should not create a systematic bias in mobility estimates because the geographic coverage of the operator is so extensive.

Quantifying Travel Using the Mobile Phone Data.

Every caller was assigned to his or her most frequently used base station/mobile phone tower on a given day, as in previous studies (30). For a given location, defined through the location of a base station, the flux of that location is defined through the number of active callers assigned to the base station and was then aggregated to the tehsil level. On an average day, ∼12% of the total population of Pakistan made a call. We measured daily travel between mobile phone towers relative to subscriber location on the previous day. Trips were aggregated to each tehsil based on the location of the origin and destination tower. We normalized trip counts by the origin tehsils’ number of active subscribers on each day (SI Text and Fig. S1). Aggregated forms of these mobility patterns can be made available upon request.
On a country-wide scale there were two significant decreases in the recorded CDR activity, one occurring during Ramadan (July 9 to August 7, 2013), reflecting less subscriber activity (Fig. S1). The decrease in the use of mobile communication services during Ramadan has been confirmed by the mobile operator as an expected effect, and this effect is seen every year. The other decrease in activity was after October 25, due to a major system upgrade to the core mobile system infrastructure, which impacted all of the customers of the mobile operator. In the collected dataset, this led to a drop in the number of location-generating subscribers recorded per day of ∼3.38 million. We adjusted the time series under the assumption that the system upgrade was a 1-d event only, based on expert opinion from the operator, and that the population in the customer base did not change its overall behavior. The flux values were adjusted after this date, assuming that the average would remain the same as in the beginning of the dataset (Fig. S1). To analyze the relationship between mobility and dengue dynamics, we approximated travel patterns between January 1 and June 1. We simulated travel, assuming that the mean number of normalized trips (normalized by flux) remained the same, and added noise [Norm(0,σ), where σ is the variance in the number of trips between all pairs of tehsils] (SI Text and Fig. S1).

Climate Data.

A. aegypti entomological and dengue viral factors are highly influenced by temperature (19). Daily mean temperature and total precipitation were recorded at 39 weather stations across Pakistan, obtained from the National Oceanic and Atmospheric Administration National Climate Data Center (Fig. S3). Temperatures peak in the middle of the year, June, July, and August; are lowest in January and December; and are variable across the country (daily averages between 9 °C and 27 °C) (for example, Fig. S3). We converted temperature to dew-point temperature based on the Numerical Terradynamic Simulation Group proposed model at the University of Montana (31). Dew-point values and temperature were then converted to relative humidity.
Fig. S3.
The temperature data from weather stations in Pakistan. (A) The location of weather stations in Pakistan. (B) Daily temperature (in Celsius) from selected weather stations in Pakistan. (C) The temperature-dependent parameters used in the ento-epidemiological framework for Karachi.

Ento-Epidemiological Framework.

To model dengue dynamics, we used an ordinary differential equation model based on a model by Lourenco and Recker to describe a dengue outbreak in Madeira, Portugal (10). This model captures the dynamics of dengue between human and mosquito hosts where the mosquito dynamics are dependent on temperature and relative humidity (SI Text and Fig. S3). Here we assume that individuals can be infected only once and we do not consider multiple serotypes. We assigned temperature-dependent epidemiological variables as in Lourenco and Recker (SI Text and Fig. S3) (10).
The relationship between relative humidity and dengue suitability has been explored in a number of environments (3, 32). We have added an additional variable to the temperature-dependent epidemiological variable, the mortality factor that is based on relative humidity (33). In contrast to ref. 10, we analyzed the reporting rate ρ along with the biting rate a (SI Text). We varied the reporting rate between 2% and 10% (a 2% reporting rate is shown) because low rates of dengue reporting have previously been found in South Asia (34). We did not fit the carrying capacity (K) explicitly and performed a sensitivity analysis, changing values of K. Although various values of K changed both a and ρ, it did not change the overall shape of the epidemic (Fig. S4 and SI Text). Using the temperature data for Karachi and the population (13.4 million), we estimated that a=0.66,ρ=0.02 (assuming K=5e6; see SI Text for the sensitivity analysis).
Fig. S4.
Sensitivity analysis on the carrying capacity (K) fitted to the biting rate (a) for each value of K, using a maximum-likelihood framework. The case data for Karachi are shown as points.

Importation of Infected Travelers from Endemic Areas in Southern Pakistan.

To estimate the number and role of importation of cases from endemic areas in southern Pakistan (Karachi) to all other tehsils, we first modeled the dengue dynamics in Karachi, using the ento-epidemiological framework described above. We are focused only on the role of travel within Pakistan as opposed to cross-border migration (Discussion). Dengue is endemic in Karachi (16), due to year-round climatic suitability and availability of vectors, and thus we assumed the date of the first introduced was day 1 (t0=1).
We then estimated the flow of infected travelers from Karachi. We first estimated the number of infected travelers per day who have left Karachi (Tt=mtβt). Based on the mobile phone data, ∼30% (in the figures 30% is shown: min, 25.6%; median, 30%; max, 34%) of subscribers have traveled outside of Karachi per day (βt) (Fig. S5). We varied this percentage between 10% and 30% to account for uncertainty in this estimate and a possible overestimate of travel because this value was based only on mobile phone subscribers (25). The daily number of infected individuals in Karachi is based on the modeled epidemic (mt). We determined the destinations of infected travelers based on the daily percentage of travelers from Karachi to all other tehsils (xt,j, mobile phone data; gt,j, gravity model to location j). We calculated these flows using either the mobile phone data or a gravity (diffusion) model.
Fig. S5.
Sensitivity analysis of estimated introductions to tehsils based on the CDR varying the percentage that travel or the probability of introduction. We analyzed the estimated introductions, using the mobile phone data to quantify travel from Karachi to other tehsils. We varied the percentage of travel from Karachi (Left) and the probability of introduction (Right).
Here, the amount of travel between two locations, i,j estimated via a basic diffusion, naive gravity, model is Ndiffusion,i,j=(popi×popj)/dist(i,j) and we estimated the population (popi,popj) of each tehsil, using Worldpop.org (www.worldpop.org), and calculated the distance (dist(i,j)) as the travel time distance between centroids of each tehsil (35). We also fitted a standard gravity model to the mobile phone data (SI Text and Fig. S2).
The raw estimates of infected travelers from Karachi to all other tehsils per day represent an upper limit. To account for variation in epidemiological and individual factors, such as the within-host viral dynamics, we assumed that a smaller subset of this upper bound of infected visitors could contribute to transmission. We sampled the actual number of effective infected travelers from a binomial distribution with a fixed probability [Pr(Ti×xt,j|γ) and Pr(Ti×gt,j|γ)]. We varied this probability between 0.001 and 0.9 (0.01 is shown); see SI Text for further details. Thus, for each day, we had an estimated number of introduced cases from Karachi to each tehsil [Pr(Ti×xt,j|γ)×mt=Yx,t, Pr(Ti×gt,j|γ)×mt=Yg,t]. The date of the first introduced case from Karachi was sensitive to both the reporting rate for Karachi and the percentage of travelers from Karachi (SI Text) (see Fig. S6 for plots of each tehsil that reported cases), although it was less sensitive to the fixed probability (Fig. S7). We also investigated the possibility of importations from Mingora into Lahore (Fig. S8) during 2013 and find that importations from Karachi to Lahore are more likely the cause of the resulting epidemic in Lahore, although these two possible causes cannot be disentangled in the case data.
Fig. S6.
Mobility models derived from mobile phone data predict the timing of introduced cases around the country that spark epidemics. Shown are the estimated introduced cases from Karachi to all tehsils that have reported at least 15 cases. The estimated introductions from the mobile phone data (blue), the diffusion model (orange), the gravity model (green), and actual case data (red) are shown. Introductions are estimated assuming that 30% of individuals travel consistent with the mobile phone data (prob = 0.01). If there were no estimated introductions from both the mobile phone data and the diffusion model, no boxes are drawn. In the majority of the KPK tehsils, both the mobile phone data and the diffusion model did not predict any introductions, likely because these cases were associated with the initial outbreak in Mingora (Alai–Peshawar). However, in the Punjab tehsils (Faisalabad, Ferozawala, Rawalpindi, and Sheikhupura) the mobile phone data consistently predicted an earlier introduction than the diffusion model and before the date of the first reported case.
Fig. S7.
The estimated introductions for Mingora/Babuzai tehsil and Lahore while varying parameters. (A and B) The estimated introductions for Mingora/Babuzai and Lahore assuming various reporting rates for Karachi. Shown are the estimated introduced cases from Karachi to Mingora and Lahore (assuming 30% of individuals travel and a probability of 0.01) with various reporting rates: 2%, 3%, 5%, and 10%. Introduced cases were estimated using the mobile phone data (blue), a gravity model fitted to the mobile phone data (orange), or a diffusion model (green). The mobile phone data consistently predict earlier introductions than either the gravity model or the diffusion model, with the diffusion model predicting the latest introductions of the three methods. In particular, the diffusion model does not report introductions into Mingora and for high reporting rates (10%), there are no predicted introductions into this tehsil. (C and D) The estimated introductions for Mingora/Babuzai and Lahore, assuming various percentages of individuals travel from Karachi. Shown are the estimated introduced cases from Karachi to Mingora and Lahore (assuming a 2% reporting rate and a probability of 0.01) with various percentages of individuals traveling from Karachi: 10%, 20%, and 30%. Unsurprisingly, as the percentage of individuals who travel increases, the number of importations to both Mingora and Lahore increases. The timing of these introductions also increases as the percentage of individuals who travel from Karachi increases. (E and F) The estimated introductions for Mingora and Lahore, assuming various probabilities. Shown are the estimated introduced cases from Karachi to Mingora and Lahore (assuming 30% of individuals travel and a 2% reporting rate) with various probabilities: 0.001, 0.01, 0.1, and 0.9 (Materials and Methods).
Fig. S8.
The estimated introductions to Lahore. Using the mobile phone data, the number of introductions from (A) Karachi and (B) Mingora to Lahore are shown. As in Fig. 3, we estimated the introduced cases from Karachi or Mingora (Swat) to Lahore. (C) A comparison of the median number of introductions per week.

Environmental Suitability and Epidemic Risk.

We defined dengue suitability as a function of temperature Zx(T) that is proportional to vectoral capacity (V) for a given location x (19) (SI Text). This measure depends on the adult vector mortality rate (determined by temperature and relative humidity) and the infectious period for adult vectors. Values of Zx(T) approximately equal to zero indicate that the environment does not permit onward transmission and can be used to determine the timing of the geographic and temporal limits of dengue transmission in Pakistan.
We created a composite measure of environmental suitability and introduction of infected travelers from Karachi, riskepidemic(x)=t=1NZx(Tt)Yx,t for a location x, where Zx(Tt) is the environmental suitability for dengue on day t. This measure sums the environmental suitability (Zx(T)) times the importation of infected travelers from Karachi per day (Yx,t).
SI Text

Mobility Dynamics Derived from Mobile Phone Data

As described in Materials and Methods, for all pairs of consecutive days over the timeframe of the data, we calculated the number of trips between all pairs of tehsils. This created a time series of travel between tehsils. We normalized this time series by the number of daily active subscribers in each origin tehsil (flux value) (Fig. S1).

Dengue in Pakistan

The first confirmed case of dengue in Pakistan occurred in 1994 (16) in Karachi. Before 2008, the majority of confirmed cases were confined to Karachi. Since then, there have been cases reported in a number of provinces with the most severe outbreak occurring in 2011 with over 20,000 confirmed cases, the majority of these in the city of Lahore in Punjab province.
The dengue season and peak have varied by year, although the peak season often occurs in the fall (October–November). A comprehensive countrywide seroprevalence survey has not been conducted in Pakistan, although smaller studies have been conducted in Karachi. In a study from 1994, DEN-1 and DEN-2 were found in Karachi (36). In a later study, authors found that DEN-2 (56.2%) was slightly more prevalent than DEN-3 (43.8%) (16) based on serum samples collected at a tertiary care hospital in Karachi from the 2006 outbreak. Based on samples collected during the 2011 outbreak in Punjab, researchers found 100% prevalence of DEN-2 in dengue-positive samples (495 total samples) (36). In the same study, 26 of these samples had concurrent infection with DEN-2 and DEN-3 with no evidence of DEN-1 or DEN-4.

Gravity Model

We also fitted a gravity model (three parameters) to the average number of trips from Karachi and all other tehsils per day (Fig. S2):
Ngm,i,j=k(popjβ)dist(i,j)γ.
Here, k=7.065,β=1.03,γ=0.156, with a reduction in deviance of 70%, using Euclidean distance (kilometers) between centroids. Using road distance (kilometers), k=7.8577,β=1.03,γ=0.149, with a reduction in deviance of 70%. Using travel time (minutes) distance, k=8.16,β=1.03,γ=0.217, with a reduction in deviance of 70%.

Converting Temperature to Relative Humidity

Relative humidity data are not regularly collected at the 39 weather stations in our dataset (see Fig. S3 for the location of weather stations). We converted temperature to relative humidity by first converting temperature to dew point (Td). We used the method first developed by the Numerical Terradynamic Simulation Group at the University of Montana (www.ntsg.umt.edu) in their MTCLIM model. Here, Td=Tnk(TdayTn), where k is a calibration coefficient (here k=1), and Tday is the estimated diurnal temperature Tday=0.45(TmaxTmean)+Tmean, where Tmax is the daily maximum temperature and Tmean is the daily mean temperature. Dew point (Td) and temperature (T) are converted to relative humidity, using RH=100exp((17.625Td)/(243.04+Td))/exp((1.7625T)/(243.04+T)).

Ento-Epidemiological Framework

To model dengue dynamics, we used an ordinary differential equation model based on a model by Lourenco and Recker to describe a dengue outbreak in Madeira, Portugal (10). This model captures the dynamics of dengue between human and mosquito hosts, where the mosquito dynamics are dependent on temperature and relative humidity (SI Text and Fig. S3). Here we assume that individuals can be infected only once and we do not consider multiple serotypes.
For a fully susceptible, fixed human population size N, when bitten by an infectious mosquito (λvh), individuals enter an incubation phase (Eh) for a mean duration of 1/γh d. Then these individuals become infectious (Ih) for 1/σh d before entering the final recovery (Rh).
The dynamics of the human population are defined by the following equations:
dShdt=λvh
dEhdt=λvhγhEh
dIhdt=γhEhσhIh
dRhdt=σhIh
N=Sh+Eh+Ih+Rh.
The vector population dynamics are originally based on a model developed by Yang et al. (37). The vectors are divided into aquatic (A, which includes eggs, larvae, and pupae) and adult females (V). Adult females are further subdivided into susceptible (SV), incubating (IV) for 1/γV d, when they become infectious (IV). The dynamics of the vector population are defined by
dAdt=θA(1AK)V(ϵA+μAV)A
dSVdt=ϵAAλhVμVVSV
dEVdt=λhVγVVEVμVVEV
dIvdt=γVEVμVVEv
V=Sv+EV+IV.
Here, ϵA denotes the rate of progression from aquatic phase to adults, μAV and μVV are the mortality rates for the aquatic phase and adults, and θA is the intrinsic oviposition rate. The additional logistic term (1A/K) is a term to measure the available capacity to receive eggs that is scaled by the carrying capacity (K).
The vector to human (λvh) and human to vector (λhv) incidence rates are assumed to be density dependent and frequency dependent, respectively,
λVh=aϕVhIVShNh
λhV=aϕhVIhSVNh,
where a is the biting rate and ϕVh,ϕhV are the vector-to-human and human-to-vector transmission probabilities per bite. The expression for dengue’s basic offspring number (Q), the mean number of viable offspring produced by one female adult during its entire lifetime, is derived as
Q=ϵAVϵAV+μAV+θVμVV.
For a given temperature T0, the expected population size of each mosquito life stage is modeled and used to initialize the system, using the temperature on the first day (T0): A(T0)=K(11/Q(T0)) and V(T0)=K(11/A(T0))ϵAV(T0)/μVV(T0).
The value of dengue’s basic reproduction number is R0=(V/Nh)a2ϕVhϕhV/μVVσh(γVV+μVV).

Temperature-Dependent Parameters

The ento-epidemiological dengue model included seven entomological parameters that were dependent on temperature (T, in degrees Celsius) (Fig. S3). These parameters included the transition rate from aquatic to adult mosquito life stages (ϵAV), mortality rate of aquatic (μAV) and adult (μVV) mosquito life stage, intrinsic oviposition rate of adult mosquito life stage (θVV), extrinsic incubation period of adult mosquito life stage (γVV), and both the human-to-vector (ϕhV) and vector-to-human (ϕVh) probabilities of transmission per infectious bite. We followed the same functional equations for these parameters as in refs. 10 and 37, which were based on temperature-controlled experiments performed on Aedes aegypti populations,
ϵAV(T)=0.1310.05723T+0.01164T20.001341T3+0.00008723T43.017×106T5+5.153×108T63.42×1010T7
μAV(T)=2.130.3787T+0.02457T26.778×104T3+6.794×106T4
μVV(T)=RHF(T,RH)(0.86920.1599T+0.01116T23.408×104T3+3.809×106T4)
θVV(T)=5.4+1.8T0.2124T2+0.01015T31.515×104T4
γVV(T)=(3.3589×103Tk)/298×exp((1,500/R)(1/2981/Tk))1+exp((6.203×1021)/R(1/(2.176×1030)1/Tk)),
where Tk is in degrees kelvin,
ϕhV(T)=1.044×103T×(T12.286)×(32.461T)1/2
ϕvh(T)=0.0729T0.97.
In addition, we have included an adult mortality factor based on relative humidity (33). Temperature and relative humidity are converted to a vapor pressure measure, VP=6.1110^(7.5T/(237.3+T)/10. This value is converted to a relative humidity factor (RHF) based on the following rules: If 10<VP<30,RHF=1.20.2VP, and if VP30,RHF=0.5.
Similar to ref. 10, we used a number of fixed parameters including the human latency period (1/γh = 2 d), the human infectious period (1/σh = 4 d), and the human population size of Lahore, Karachi, or Mingora. We performed a sensitivity analysis on the mosquito carrying capacity (K).
We fitted the biting rate and the reporting rate, using a maximum-likelihood framework.

Fitting the First Introduced Case into Lahore and Mingora from Karachi, Using the Epidemiological Data

To compare our importation estimates from the mobile phone data, we estimated the first introduced case into both Mingora (Babuzai tehsil) and Lahore, using only the dengue case data for these locations. Using the same ento-epidemiological framework, we estimated the date of the first introduced case to Mingora or Lahore (T0) for a range of reporting rates (ρ). We used the same biting rate and carrying capacity for Karachi, although we recalculated the seven entomological parameters based on the daily temperatures in Mingora or Lahore. For both locations, we varied the reporting rate (ρ=0.03,0.06,0.1,0.3,0.5,0.6,0.7,0.8,0.9) and estimated T0 using a maximum-likelihood framework. In Mingora, we estimate T0 is between days 210 and 231. For comparison, day 231 (August 19) is the date of the first reported case in Mingora. In Lahore, we estimate T0 is between days 124 and 242. For comparison, 95% of the cases occurred after day 259 (September 16), which is ∼4 mo after day 124 and ∼2 wk after day 242.
Using a range of reporting rates (for comparison the fitted reporting rate for Karachi is 0.03), we estimated the date of the first introduced case (T0), using the ento-epidemiological framework described in Materials and Methods. For low reporting rates (below 30%), T0 increases with the reporting rate. For high reporting rates, the values remain constant at 231 (Mingora) and 130 (Lahore).
Reporting rate, ρMingora, estimated T0Lahore, estimated T0
0.03210124
0.06221130
0.1228130
0.3231130
0.5231130
0.6231130
0.7231130
0.8231130
0.9231130

Introductions into Other Tehsils from Karachi

We simulated the introduction of infected travelers from Karachi to all other tehsils based on either the mobile phone data or a gravity model. We first estimated the number of infected individuals in Karachi each day, using the ento-epidemiological framework described above. In each simulation we selected the percentage of infected individuals in Karachi who traveled to another tehsil (range of 10%, 20%, or 30%). This percentage provided an estimate on the number of infected travelers from Karachi. The destination for each infected traveler was determined based on the actual distribution of trip destinations from Karachi per day. We then varied the probability that these infected visitors in all other tehsils would impact the dengue dynamics in that location and sampled a binomial distribution with a fixed probability (probability varied from 0.001 to 0.9). For each fixed percentage of infected individuals who travel and probability (P), we ran 200 simulations.

Introductions from Swat to Lahore

The Swat epidemic peaked 2 mo earlier than the Lahore epidemic; we also investigated the possible role of importation from Swat to Lahore. First we compared travel from Swat tehsils that have had at least 15 cases of dengue (Alai, Alpuri, Babuzai, Bahrain, Barikot, Charbagh, Chiksar, Dassu, Kabal, Kalam, Khwazakhela, Mamund, Matta Khararai, Matta Sebujni, and Peshawar) to the amount of travel from Karachi to Lahore. Based on the CDR data, the amount of travel from Karachi to Lahore is an order of magnitude greater than the amount of travel from Swat to Lahore (number of trips per day from Karachi to Lahore, median, 18,980; min, 10,130; max, 21,940) (number of trips per day from Swat to Lahore, median, 4,463; min, 2,556; max, 5,345).
We modeled the dengue dynamics in Swat, using the same ento-epidemiological framework described in Materials and Methods with ρ=0.03, T0=210,a=0.66.

A. aegypti Suitability in Pakistan

A. aegypti have been observed in many parts of Pakistan, including at least 13 major cities, based on a genetic population structure study conducted in 2009 and 2010 (38). In Karachi, the vector has been present in all parts of the city since at least 2000 (39), with a follow-up study in 2010 (40). It has been hypothesized that malaria control programs in the 1950s reduced the A. aegypti population and since the early 2000s the tire trade has spread the vector to other parts of Pakistan from Karachi.

Dengue Temperature Suitability

We used a relative dengue suitability index, a function proportional to vectoral capacity. This measure is based on temperature and relative humidity and does not take into account other factors such as population density or percentage urban, which has been shown to impact Aedes environmental suitability (19).
We defined dengue suitability as a function of temperature Zx(T) that is proportional to vectoral capacity (V) for a given location x (19). Values of Zx(T) approximately equal to zero indicate that the environment does not permit onward transmission and can be used to determine the timing of the geographic and temporal limits of dengue transmission in Pakistan,
V=ma2μVeμVVγVV,
where m is the ratio of adult female mosquitoes to humans, a is the human blood feeding rate, μVV is the instantaneous per-capita death rate of adult mosquito females, and γVV is the length of the extrinsic incubation period for the pathogen (complete equations in SI Text). In our model, μVV is also scaled by a mortality factor based on absolute humidity. Letting λ=mμVV, Zx(T)=expμVV(T)γVV(T)/μVV(T)2V(T)=λa2(eμVV(T)γVV(T)/μVV(T)2), where Z(T) effectively measures the relative number of infectious mosquitoes supported in an environment with temperature T, given a constant λ,a. Given the poor epidemiological data to estimate λ,a in many locations in Pakistan, Zx(T) provides a relative measure of vectoral capacity for a location x and temperature T.

Acknowledgments

A.W. is funded by a James S. McDonnell Foundation postdoctoral fellowship. M.F.B. is supported by Wellcome Trust Grant 098511/Z/12/Z. C.O.B., A.W., and M.A.J. were supported by the Models of Infectious Disease Agent Study program (Cooperative Agreement 1U54GM088558).

Supporting Information

Supporting Information (PDF)
Supporting Information

References

1
; WHO Global Strategy for Dengue Prevention and Control (World Health Organization, Geneva, 2012).
2
MG Guzman, E Harris, Dengue. Lancet 385, 453–465 (2015).
3
S Bhatt, et al., The global distribution and burden of dengue. Nature 496, 504–507 (2013).
4
C Heintze, M Velasco Garrido, A Kroeger, What do community-based dengue control programmes achieve? A systematic review of published evaluations. Trans R Soc Trop Med Hyg 101, 317–325 (2007).
5
ME Halloran, Jr IM Longini, Emerging, evolving, and established infectious diseases and interventions. Science 345, 1292–1294 (2014).
6
KE Jones, et al., Global trends in emerging infectious diseases. Nature 451, 990–993 (2008).
7
ST Stoddard, et al., House-to-house human movement drives dengue virus transmission. Proc Natl Acad Sci USA 110, 994–999 (2013).
8
MR Nunes, et al., Air travel is associated with intracontinental spread of dengue virus serotypes 1-3 in Brazil. PLoS Negl Trop Dis 8, e2769 (2014).
9
CP Simmons, JJ Farrar, V Nguyen, B Wills, Dengue. N Engl J Med 366, 1423–1432 (2012).
10
J Lourenço, M Recker, The 2012 Madeira dengue outbreak: Epidemiological determinants and future epidemic potential. PLoS Negl Trop Dis 8, e3083 (2014).
11
ST Stoddard, et al., The role of human movement in the transmission of vector-borne pathogens. PLoS Negl Trop Dis 3, e481 (2009).
12
GMSS Vazquez-Prokopec, et al., Usefulness of commercially available GPS data-loggers for tracking human movement and exposure to dengue virus. Int J Health Geogr 8, 68 (2009).
13
DM Morens, AS Fauci, Emerging infectious diseases: Threats to human health and global stability. PLoS Pathog 9, e1003467 (2013).
14
YL Hii, H Zhu, N Ng, LC Ng, J Rocklov, Forecast of dengue incidence using temperature and rainfall. PLoS Negl Trop Dis 6, e1908 (2012).
15
A Amarasinghe, GW Letson, Dengue in the Middle East: A neglected, emerging disease of importance. Trans R Soc Trop Med Hyg 106, 1–2 (2012).
16
SB Rasheed, RK Butlin, M Boots, A review of dengue as an emerging disease in Pakistan. Public Health 127, 11–17 (2013).
17
A Ali, et al., Seroepidemiology of dengue fever in Khyber Pakhtunkhawa, Pakistan. Int J Infect Dis 17(7):e518–523. (2013).
18
N Bharti, et al., Explaining seasonal fluctuations of measles in Niger using nighttime lights imagery. Science 334, 1424–1427 (2011).
19
OJ Brady, et al., Global temperature constraints on Aedes aegypti and Ae. albopictus persistence and competence for dengue virus transmission. Parasit Vectors 7, 338 (2014).
20
J Aldstadt, et al., Space-time analysis of hospitalised dengue patients in rural Thailand reveals important temporal intervals in the pattern of dengue virus transmission. Trop Med Int Health 17, 1076–1085 (2012).
21
; United Nations, Department of Economic and Social Affairs, Population Division, Trends in International Migrant Stock: The 2013 Revision—Migrants by Destination and Origin (United Nations database, POP/DB/MIG/Stock/Rev.2013/Origin). Available at www.un.org/en/development/desa/population/publications/pdf/migration/migrant-stock-origin-2013.pdf. Accessed June 1, 2015. (2013).
22
N Gupta, S Srivastava, A Jain, UC Chaturvedi, Dengue in India. Indian J Med Res 136, 373–390 (2012).
23
S Lai, et al., The changing epidemiology of dengue in China, 1990-2014: A descriptive analysis of 25 years of nationwide surveillance data. BMC Med 13, 100 (2015).
24
M Mardani, F Abbasi, M Aghahasani, B Ghavam, First Iranian imported case of dengue. Int J Prev Med 4, 1075–1077 (2013).
25
A Wesolowski, N Eagle, AM Noor, RW Snow, CO Buckee, The impact of biases in mobile phone ownership on estimates of human mobility. J R Soc Interface 10(81):20120986. (2013).
26
N Imai, I Dorigatti, S Cauchemez, NM Ferguson, Estimating dengue transmission intensity from sero-prevalence surveys in multiple countries. PLoS Negl Trop Dis 9, e0003719 (2015).
27
PS Wikramaratna, CP Simmons, S Gupta, M Recker, The effects of tertiary and quaternary infections on the epidemiology of dengue. PLoS One 5, e12347 (2010).
28
NS Reich, et al., Interactions between serotypes and dengue highlight epidemiological impact of cross-immunity. J R Soc Interface 10(86):20130414. (2013).
29
CN Koo, et al., Evolution and heterogeneity of multiple serotypes of dengue virus in Pakistan, 2006-2011. Virol J 10(1):275. (2013).
30
A Wesolowski, et al., Quantifying the impact of human mobility on malaria. Science 338, 267–270 (2012).
31
E Eccel, Estimating air humidity from temperature and precipitation measures for modelling applications. Meteorol Appl 19, 118–128 (2012).
32
Hales S, deWet N, Maindonald J, Woodward A (2002) Potential effect of population and climate changes on global distribution of dengue fever: An empirical model. Lancet 360(9336):830–834.
33
DA Focks, DG Haile, E Daniels, GA Mount, Dynamic life table model for Aedes aegypti (Diptera: Culicidae): Analysis of the literature and model development. J Med Entomol 30, 1003–1017 (1993).
34
Organization WH (2010) Situation Update of Dengue in the SEA Region 2010 (World Health Organization, Geneva).
35
A Wesolowski, et al., The use of census migration data to approximate human movement patterns across temporal scales. PLoS One 8, e52971 (2013).
36
DSIA Akram, A Igarashi, T Takasu, Dengue virus infection among children with undifferentiated fever in Karachi. Indian J Pediatr 65, 735–740 (1998).
37
HM Yang, ML Macoris, KC Galvani, MT Andrighetti, DM Wanderley, Assessing the effects of temperature on the population of Aedes aegypti, the vector of dengue. Epidemiol Infect 137, 1188–1202 (2009).
38
SB Rasheed, M Boots, AC Frantz, RK Butlin, Population structure of the mosquito Aedes aegypti (Stegomyia aegypti) in Pakistan. Med Vet Entomol 27, 430–440 (2013).
39
RMZS Tariq, Why the population of dengue vector mosquitoes is increasing day-by-day in Karachi and other areas of Sindh, Pakistan? Pak J Entomol 15, 7–10 (2000).
40
RMAI Tariq, SS Qadri, Population dynamics and mechanical control of dengue vector mosquitoes Aedes aegypti and Aedes unilineatus in seven towns of Karachi. Pak J Entomol 25, 21–26 (2010).

Information & Authors

Information

Published in

Go to Proceedings of the National Academy of Sciences
Go to Proceedings of the National Academy of Sciences
Proceedings of the National Academy of Sciences
Vol. 112 | No. 38
September 22, 2015
PubMed: 26351662

Classifications

Submission history

Published online: September 8, 2015
Published in issue: September 22, 2015

Keywords

  1. dengue
  2. human mobility
  3. Pakistan
  4. mobile phones
  5. epidemiology

Acknowledgments

A.W. is funded by a James S. McDonnell Foundation postdoctoral fellowship. M.F.B. is supported by Wellcome Trust Grant 098511/Z/12/Z. C.O.B., A.W., and M.A.J. were supported by the Models of Infectious Disease Agent Study program (Cooperative Agreement 1U54GM088558).

Notes

This article is a PNAS Direct Submission.

Authors

Affiliations

Amy Wesolowski
Department of Epidemiology, Harvard T. H. Chan School of Public Health, Boston, MA 02115;
Center for Communicable Disease Dynamics, Harvard T. H. Chan School of Public Health, Boston, MA 02115;
Taimur Qureshi
Telenor Research, Telenor Group, N-1360 Fornebu, Norway;
Maciej F. Boni
Oxford University Clinical Research Unit, Wellcome Trust Major Overseas Programme, Ho Chi Minh City, Vietnam;
Centre for Tropical Medicine, Nuffield Department of Clinical Medicine, University of Oxford, Oxford OX3 7FZ, United Kingdom;
Pål Roe Sundsøy
Telenor Research, Telenor Group, N-1360 Fornebu, Norway;
Michael A. Johansson
Center for Communicable Disease Dynamics, Harvard T. H. Chan School of Public Health, Boston, MA 02115;
Division of Vector-Borne Diseases, Centers for Disease Control, San Juan, Puerto Rico 00920;
Syed Basit Rasheed
Department of Zoology, University of Peshawar, Peshawar 25120, Pakistan
Kenth Engø-Monsen
Telenor Research, Telenor Group, N-1360 Fornebu, Norway;
Caroline O. Buckee1 [email protected]
Department of Epidemiology, Harvard T. H. Chan School of Public Health, Boston, MA 02115;
Center for Communicable Disease Dynamics, Harvard T. H. Chan School of Public Health, Boston, MA 02115;

Notes

1
To whom correspondence should be addressed. Email: [email protected].
Author contributions: A.W. and C.O.B. designed research; A.W., T.Q., M.F.B., K.E.-M., and C.O.B. performed research; A.W., M.F.B., M.A.J., S.B.R., and K.E.-M. contributed new reagents/analytic tools; A.W., T.Q., P.R.S., S.B.R., and K.E.-M. analyzed data; and A.W., T.Q., M.F.B., P.R.S., M.A.J., S.B.R., K.E.-M., and C.O.B. wrote the paper.

Competing Interests

Conflict of interest statement: M.F.B. has worked as a paid consultant to Visterra, Inc. in Cambridge, MA.

Metrics & Citations

Metrics

Note: The article usage is presented with a three- to four-day delay and will update daily once available. Due to ths delay, usage data will not appear immediately following publication. Citation information is sourced from Crossref Cited-by service.


Citation statements




Altmetrics

Citations

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

Cited by

    Loading...

    View Options

    View options

    PDF format

    Download this article as a PDF file

    DOWNLOAD PDF

    Get Access

    Login options

    Check if you have access through your login credentials or your institution to get full access on this article.

    Personal login Institutional Login

    Recommend to a librarian

    Recommend PNAS to a Librarian

    Purchase options

    Purchase this article to get full access to it.

    Single Article Purchase

    Impact of human mobility on the emergence of dengue epidemics in Pakistan
    Proceedings of the National Academy of Sciences
    • Vol. 112
    • No. 38
    • pp. 11739-E5378

    Media

    Figures

    Tables

    Other

    Share

    Share

    Share article link

    Share on social media