Ethics statement

All survey respondents provided verbal informed consent. The project was approved by the Institutional Review Board of NIPSOM (Memo number: NIPSOM/IRB/2017/320/1) and the Bangladesh Medical Research Council National Research Ethics Committee (Registration number: 06310082017). The secondary data analysis of the de-identified survey data was determined by the Harvard T H Chan School of Public Health Institutional Review Board (Protocol number: IRB18-1693) as exempt.

Epidemiological data

A cluster-randomized household survey was conducted in Dhaka city in July 2017 during the peak of the chikungunya outbreak. One hundred clusters were defined in Dhaka city based on administrative boundaries. Approximately 30 households were randomly selected within each cluster, resulting in a total of 3253 households. One individual within each household was asked about symptoms of chikungunya (fever, joint pain, and rash), and the date of onset of the symptoms. Of the suspected cases that met the clinical definition (2518 individuals), a subset were randomly selected for laboratory confirmation of infection with the chikungunya virus (CHIKV). Serum samples were obtained from 1487 individuals and an immunochromatographic test was used to check for the presence of antibodies—immunoglobulin M (IgM) and immunoglobulin G (IgG)—against CHIKV for a random subset of samples (1286 samples). IgM antibodies for CHIKV can be detected 3 to 4 days after clinical onset of illness, and remain detectable for 1 to 3 months post-infection [17, 28]. IgG antibodies, on the other hand, can remain detectable for years [17, 29]. We define an individual to have a confirmed recent infection in 2017, if their sample tests positive for IgM antibodies, as it is unlikely they were infected in a previous outbreak. The delay between infection and the ability to detect IgM antibodies could lead to an under estimation of cases. Due to the retrospective nature of the survey, and the fact that the majority of those reporting symptoms did so more than a week before the survey was conducted, we expect the under-estimation due to the delay to be minimal (and only affecting those who reported symptoms in the final week of the survey). Of the 1286 samples tested for IgM, 895 were positive. Concurrent to the household survey, an entomological team also collected larvae and mosquito samples from pots and containers inside, and in the vicinity of, the house.

Mobility data

We derived mobility estimates from CDR data from Telenor Group’s mobile operator Grameenphone in Bangladesh, with over 64 million subscribers, using methods that have been previously described [3, 4, 8]. The anonymized mobile phone data, covering the the time period between April 1, 2017 and September 30, 2017, was used to estimate movement using the location information encoded in the CDR. Subscribers were assigned to a tower location for each day in the dataset according to their most frequently used routing tower for calls. Trips were calculated based on changes in a subscriber’s assigned tower location from the previous day. Due to data regulations and privacy concerns, all data were aggregated temporally (daily) and spatially (to the smallest administrative level, a union) based on tower locations, as described previously [4, 8, 30]. The aggregated data thus consisted of the daily number of subscribers for each union and the total number of trips between all pairs of unions for a 6 month period.

Modeling drivers of spatial heterogeneity within Dhaka

To understand the epidemiology of the outbreak within Dhaka, we first examined the factors associated with spatial heterogeneities in chikungunya incidence within Dhaka. We aggregated the survey data to 95 administrative units for which we also had information on socioeconomic and environmental indicators (S1 Fig). This spatial aggregation was necessary for correctly matching our data sources, as union boundary definitions have changed several times in Dhaka. To account for the fact that different proportions of suspected cases were sent for lab confirmations in different locations, we estimated the number of chikungunya cases in each location, O_i, based on the number of confirmed cases, C_i, and suspected cases (reported symptoms but not tested) S_i, in each location, and the global ratio of confirmed cases to the number tested, T_i:

We used a generalized linear model to assess the impact of demographic, socioeconomic, and environmental factors on estimated cases in each survey location. Specifically, we assumed that our outcome of interest, O_i follows a Poisson distribution with mean, η_i:

We modeled the logarithm of η_i as a linear combination of a set of independent variables, with the number of households, households_i as an offset. The full model specification is:

Human-mosquito transmission model

Since CHIKV is known to cause asymptomatic infections [17], the survey data cannot be used directly to obtain the total burden of chikungunya in Dhaka, since only symptomatic individuals were tested for infection. This is particularly important for estimating importation of cases from Dhaka, as asymptomatic individuals are more likely to travel than symptomatic individuals given the nature and severity of symptoms. We used a mechanistic framework (S4 Fig) to model chikungunya dynamics in Dhaka to estimate daily prevalence of symptomatic and asymptomatic infections. Due to the model complexity, we estimate a single model for all of Dhaka city, rather than a spatially-explicit model.

In this compartmental framework, susceptible humans are exposed to the pathogen through the bite of an infected mosquito. After a latent period, people are either asymptomatically infectious or symptomatically infectious, during which time they can infect susceptible mosquitoes, before recovering and acquiring life-long immunity. In this model, S is the proportion of susceptible humans, E is the proportion of humans in the latent incubation period, I is the proportion that are symptomatically infections, I^A is the proportion that are asymptomatically infectious, and R is the proportion that have recovered. Since infection with CHIKV is thought to confer lifelong immunity [17], recovered individuals remain in the R compartment. For the purpose of fitting to the data, we also tracked the cumulative proportion of symptomatic, C, and asymptomatic, C^A, infections. In the mosquito population, S^M is the proportion susceptible, E^M is the proportion in the latent phase and, I^M is the proportion infected. The human-mosquito dynamics are modeled as:

Table 1 describes the model parameters.

Table 1

Parameter definitions for the chikungunya transmission model.

Parameter	Definition	Value	Reference
β1	mosquito-to-human transmission rate	estimated
ϕ	proportion of infected people who become symptomatic	estimated
β2	human-to-mosquito transmission rate	estimated
1/λ1	average latent period of infection in humans	4 days	[17, 34–38]
1/γ	average duration of infectiousness	4 days	[17, 35, 36, 38]
1/μ	average mosquito lifespan	15 days	[35, 36]
1/λ2	average latent period of infection in mosquitoes	3 days	[34–36]

The basic reproduction number, R₀, is the spectral radius of the next-generation matrix [36]:

Note, that the expected number of secondary infections in humans that result from a single infected human, ($$ R_0^H$$) is $$ R_0^2$$, since two generations are required to transmit an infection from human to human [39]. Since we are primarily interested in human to human infection, we report our estimate of $$ R_0^H$$.

Based on the range of estimates from the literature (see Table 1), we fixed the average latent period of infection in humans, 1/λ₁, at 4 days, average latent period of infection in mosquitoes, 1/λ₂, at 3 days, and the average duration of infectiousness in humans, 1/γ, at 5 days. We assumed a stable mosquito population, such that the mosquito birth and death rates, μ, are the same, and equal to the inverse of the mosquito lifespan, which is fixed at 15 days. We assumed that 80% of the population were susceptible at the start of the outbreak, to account for the fact that Bangladesh had previously experienced chikungunya outbreaks in 2008, 2011 and 2012 [21–23], although these outbreaks were primarily reported in a few small rural communities outside Dhaka. In the survey population, about 5% of individuals had been infected with CHIKV prior to 2017 (IgG positive and IgM negative), although the actual proportion is likely to be higher since only individuals who reported having recent symptoms were tested. Here, we assume that 20% of the population had been infected previously, but varied this proportion in additional sensitivity analyses. The first suspected case in the survey occurred on April 2nd, and the virus was likely circulating prior to that. Thus, the initial incidence in the model was fixed according to the first observation in the survey data. The mosquito-to-human transmission rate, β₁, the human-to-mosquito transmission rate, β₂, and the proportion of infected people who become symptomatic, ϕ, were estimated by fitting the simulated incidence to the survey data.

The weekly incidence in the survey data is given by, $$ I_w^{obs}=O_w/population$$, where O_w was calculated using Eq 1 for each week, w, of the outbreak, and population is the total number of individuals surveyed, adjusting for the fact that households were not all surveyed simultaneously. To compare the model output to the survey data, we tracked the daily incidence of symptomatic infections ($$ I_t^{sim}=C_{t+1}-C_t$$). We assumed that there is a fixed lag of 2 days for infected individuals to develop symptoms [36]. Thus, the daily incidence of symptomatic infections was assumed to be equivalent to $$ I_t^{sim}$$ from 2 days ago. The daily symptomatic incidence was summed every 7 days to obtain weekly incidence, $$ I_w^{sim}$$. The maximum-likelihood parameter values were estimated by assuming normally-distributed residuals between the observed data and model predictions [36, 38]. We also fit the model assuming a measurement error structure instead i.e. we maximize the likelihood that $$ I_w^{sim}$$ is drawn from a Poisson distribution with mean given by the observed data, and obtained qualitatively similar results (S5 Fig). The log-likelihood was maximized through the Nelder-Mead algorithm within the optim package in R. We repeated the optimization 100 times with starting values for each parameter drawn from reasonable parameter ranges (0.2 ≤ β₁ ≤ 0.8;0.5 ≤ ϕ ≤ 1;0.2 ≤ β₂ ≤ 0.8), using latin hypercube sampling to explore the parameter space. The best-fit model was the set of parameter estimates that had the highest log-likelihood value. We show, in additional parameter identifiability analysis with simulated data, that joint estimation of these parameters is possible (S6 Fig) with the observed data.

We calculated uncertainty bounds for parameter estimates by assuming an observational error structure for the incidence timeseries. Specifically, we calculated approximate confidence intervals for the parameter estimates by simulating 200 realizations of the best-fit weekly incidence curve using parametric bootstrap with a Poisson error structure, as has been described in previous studies [40–42]. The best fit values of the weekly incidence, $$ I_t^{sim}$$, was perturbed by adding a simulated error. For each value of $$ I_t^{sim}$$, we simulated $$ I_t^{\prime sim}$$ which was drawn from a Poisson distribution with mean, $$ I_t^{sim}$$. Parameters were estimated for each of the 200 simulated realizations of weekly incidence. The nominal 95% confidence intervals were constructed from the distribution of the parameter estimates. This approach for constructing the confidence intervals assumes that the only source of noise is observational error.

To test model sensitivity to parameter values, we conducted additional sensitivity analyses using latin hypercube sampling to explore the parameter space and computed partial rank correlations between model parameters and model output. The model output was sensitive to different parameter values (S7A and S7C Fig). However, the partial rank correlation coefficients suggest that the model output variations were largely driven by changes in the values of the parameters we estimated (β₁, β₂, and ϕ), rather than the fixed parameters (S7B and S7D Fig). The maximum-likelihood parameter estimates were stable across a range of initial starting conditions, and the confidence intervals, constructed assuming an observational error structure, were relatively tight. Further, the parameter estimates were stable across a range of values for the initial proportion susceptible (S1 Table), the only fixed parameter for which we had no information from previous studies.

Quantifying importation of infected travelers from Dhaka

We used the modeled dynamics of chikungunya in Dhaka, and the travel pattern derived from the mobile phone data to estimate the daily number of infected travelers from Dhaka to the rest of Bangladesh, using methods described in [8]. We first estimated the daily number of infected travelers leaving Dhaka, T_t, as T_t = m_t × π_t, where m_t is the product of population size and prevalence from the best-fit model (described above) and π_t is the daily proportion of people traveling out of Dhaka. To quantify where people are traveling to, we calculated the daily proportion of travelers from Dhaka who travel to each location, x_t,j.

We quantified importation using mobility estimates from the mobile phone data, as well as a diffusion (gravity) model for comparison. For the model parameterized by the mobile phone data, we estimated π_t as the daily proportion of mobile phone subscribers who travel out of Dhaka, and we assume x_t,j is equal to the proportion of subscribers leaving Dhaka who go to location, j. For the diffusion model, we use a gravity model parameterized with estimates from the literature [43]. We calculated the number of travelers to each location as, ($$ X_{t,j}=k(p^{\alpha }\times p_j^{\beta })/d_j^{\gamma }$$; where p is the population in Dhaka, p_j is the population in location j, d_j is the Euclidean distance from Dhaka to location j, and α, β, and γ are parameter estimates from [43]. In the diffusion model, π_t was given by ∑_∀j X_t,j/p. The proportion traveling to each location, x_t,j, was given by X_t,j/∑_∀j X_t,j.

Since not all importations will lead to local transmission, the number of effective infected travelers to each location was sampled from a binomial distribution, similar to a previous study on dengue importation [3], where the probability of “success” (i.e. having an effective infected traveler), p_eff, was varied from 0.01 to 0.5, and the number of trials was equal to T_t × x_t,j. Values of p_eff below one reflect the fact that heterogeneities in epidemiological, environmental, and individual factors are likely to result in a smaller number of effective infected travelers than T_t × x_t,j. In the main results, we assume that on average, 10% of importations result in local transmission, but vary this percentage in sensitivity analyses (S8–S10 Figs).

Epidemiology of chikungunya in Dhaka

We estimated a total of 1589 cases of symptomatic chikungunya infection in the 3253 surveyed individuals, based on the proportion of tested cases that were positive for chikungunya and the total number of individuals that reported symptoms (Fig 1A). There was significant spatial heterogeneity in incidence of symptomatic cases within Dhaka. The estimated incidence ranged from 46 to 824 cases per 1000 people across the 95 administrative units surveyed (Fig 1B). Our results indicate that both the abundance of the vector and locals demographics were important drivers of the dynamics of the outbreak. As expected, the Breteau Index, a measure of vector abundance, showed substantial spatial variation (S11 Fig), and had the largest impact—a one unit change in the BI was associated with a 9% [95% CI: 2.5%, 15.8%] increase in symptomatic incidence (Fig 1C). We also found statistically significant positive associations between the incidence rate and population density and connectedness (as measured by the mobility network degree).

Fig 1

(A) Weekly number of confirmed and estimated cases in the survey data. Orange dashed line shows the number of cases reported to the Institute of Epidemiology, Disease Control and Research (IEDCR) over the course of the outbreak. (B) Estimated cases per household for each survey location. (C) Estimated relative risks (incidence rate ratios) from poisson regression model. 95% Confidence intervals were calculated by applying the delta method [44] to robust standard errors.

We estimated the daily incidence of symptomatic and asymptomatic cases in Dhaka over the course of the outbreak by fitting a mechanistic human-mosquito transmission model to the observed incidence time series. We estimated the following model parameters by fitting the simulated incidence to the survey data: the mosquito-to-human transmission rate, β₁, the human-to-mosquito transmission rate, β₂, and the proportion of infected people who become symptomatic, ϕ. Maximum-likelihood parameter estimates, and associated uncertainties, are shown in Table 2.

Table 2

Parameters were estimated by fitting the mechanistic chikungunya model to observed incidence from the survey.

Parameter estimates for the best-fit model.

Parameter	Maximum-likelihood estimate	Approximate 95% Cl
β1	0.51	(0.32, 0.75)
ϕ	0.78	(0.73,0.85)
β2	0.16	(0.10,0.29)
$$	4.20	(3.83, 4.62)

The observed survey incidence had two peaks, which may have been driven by the timing of the monsoon rainfall (S12 Fig) [45] or may be due to patchy reporting of cases. Reproducing the bi-modal nature of the outbreak, however, is not possible without increasing model complexity and adding external forcing in the model. Nonetheless, our best-fit model accurately reproduces the general trend in the observed incidence and the size of the outbreak (Fig 2). The simulated attack rate for symptomatic cases, using the best-fit model parameters, from the start of the outbreak to the end of July, was 51% (standard deviation across simulations with observational error: 1.8), which compares well with the 49% estimated from the survey. According to the best-fit model, 78% [95% CI: 73, 85] of the infections were symptomatic. We estimated a human-mosquito transmission rate of 0.51 [95% CI: 0.32, 0.75]; this is equivalent to the number of mosquito bites per human per day resulting in transmission (accounting for imperfect transmission of the pathogen). Our estimates for both β₁ and β₂ are similar to previous estimates for chikungunya using a similar model [36].

Fig 2

(A) Observed and fitted weekly incidence. The blue line shows the observed incidence from the survey. The orange dotted line shows the best-fit model. The shaded region shows the full range of incidence for 100 simulations assuming a poisson error structure. (B) Daily prevalence of symptomatic and asymptomatic infections simulated using the best-fit model parameters. The grey shaded region indicates the time period beyond the survey timeframe.

Using the best-fit model parameters, we estimated the human-to-human basic reproduction number, $$ R_0^H$$, to be 4.20 [95% CI: 3.83, 4.62], which is similar to previous estimates of R₀ for chikungunya [35, 46, 47]. We estimated a peak prevalence of 47 cases per 1000, with the peak incidence occurring in the first week of July. Interestingly, these estimates are far higher than the officially reported number of 13,176 total cases between April and September [45].

Population mobility in Bangladesh and predicted importation from Dhaka

Fig 3 characterizes the population movement patterns in Dhaka, and travel from Dhaka to the rest of the country, derived from geo-located, anonymized and aggregated CDR data. We restricted the analysis to the area covered by the epidemiological survey, namely the core of Dhaka city. The mobile phone data for Dhaka represents about 4.4 million subscribers within the survey area, compared to a population size of 7.3 million. Within Dhaka, we found spatial heterogeneity in both mobility flux, a measure of average daily flow in and out of a location, and mobility network degree, a measure of the connectedness of each union. Specifically, neighborhoods in the northwest and southwest were more connected and had larger flows in and out of the area, proportional to their population size, compared to other locations. Increased connectivity, specifically having an additional degree in the mobility network, was associated with a 7% [95% CI: 0.5%, 14.5%] increase in the incidence rate of symptomatic cases (Fig 1C).

Fig 3

Population movement patterns in Dhaka.

(A) Number of subscribers in a union versus the population size, for the unions covered by the survey. Map of Dhaka showing two measures of mobility: (B) the number of degrees in the network that had above average weight and (C) the average flow of people in and out of an area as a proportion of the population size. (D) Top panel: The daily proportion of subscribers who traveled from Dhaka to the rest of the country. This data represents the 4.4 million subscribers within the study area. Bottom panel: Daily number of subscribers (log count) traveling from Dhaka to all other districts. Districts are arranged from top to bottom by distance from Dhaka in ascending order.

Large volumes of travel occurred daily between Dhaka and all other regions in the country, particularly on weekends and holidays. On average, 6% of subscribers within our study area traveled to other parts of the country daily. Fig 3D shows the daily proportion of subscribers leaving Dhaka. Weekends (Fridays and Saturdays) and the Eid holidays (which take place twice a year) can be clearly identified. During the start of the first Eid holidays, on June 23rd, approximately 1.7 million subscribers (40% of total subscribers) left Dhaka to visit other parts of the country.

We combined our mobility and prevalence estimates to predict the importation of chikungunya from Dhaka to the rest of the country. We compare our estimates to a standard diffusion model, to provide a comparison with predictions that would be made without mobile phone data, and where possible we compare our predictions to places that did report CHIKV cases during the outbreak, though these were limited. Our results highlight important deviations in actual travel patterns compared to predictions from a diffusion model. We predicted a high likelihood of importation into areas with large volumes of travel from Dhaka, particularly during the Eid holidays. Compared to the naive diffusion model, the model parameterized with mobile phone data showed much greater variation in imported cases across space and time (Figs 4 and 5). Specifically, for some districts that reported suspected cases, such as Bogra and Dinajpur, we predicted a larger number of imported cases than would be expected based on its distance from Dhaka. In the diffusion model importation decayed rapidly with distance, with very few imported cases to locations far from Dhaka. In contrast, the mobility estimates showed high volumes of travel to some locations, such as Dinajpur district, that are far from Dhaka. In general, the diffusion model predicted later introduction to places far from Dhaka, and earlier introduction to places closer to Dhaka, compared to the model with mobile phone data (Fig 4).

Fig 4

Spatial heterogeneity in estimated importation time.

Introduction is defined here as the importation of at least ten cases, and measured as days since the start of the outbreak in Dhaka. The map on the left shows the estimated introduction time for each district based on mobility estimates. Lighter colors indicate earlier introduction. The districts that reported suspected cases to IEDCR are highlighted in red. The location of Dhaka city is indicated with the black circle. The difference between estimates using mobile phone data and a diffusion model are highlighted in the map on the right. In general, the diffusion model predicts earlier introduction (green) to nearby locations, and later (purple) or no introduction to more distant locations.

Fig 5

Daily simulated imported cases (log count) from Dhaka to other districts in Bangladesh using (A) mobility estimates versus a (B) diffusion model.

Importations are simulated assuming only the asymptomatic infected are traveling (p_eff = 0.1). Each row represents a district; districts are arranged from top to bottom by distance from Dhaka in ascending order. The districts that reported suspected cases to IEDCR are highlighted in red.

The Eid holidays were associated with a substantial increase in the estimated number of imported cases, which was not captured in the diffusion model (Fig 5). This reflects the large volume of travel from Dhaka to the rest of the country during the start of the Eid holidays. Our estimates of the number of imported cases is conservative, as we are only considering travel from the 95 unions included in the survey. In reality, the outbreak in Dhaka likely affected all unions in the greater metropolitan area, thus creating a much larger pool of potential infected travelers, especially during Eid.

Due to the lack of a surveillance and disease reporting system, only a handful of districts outside Dhaka reported cases to the Institute of Epidemiology, Disease Control and Research (IEDCR). Since reporting was not mandatory, there is no information about chikungunya incidence in the majority of districts in Bangladesh. However, the geographic distribution of the reporting districts (Fig 4 and S13 Fig) suggests that the outbreak had spread to many parts of the country, including districts far from Dhaka such as Dinajpur.

Eight of the 64 districts in Bangladesh, reported at least one probable chikungunya case, and 17 reported suspected cases, with the majority of cases reported in Narshingdi, a district close to Dhaka where we predicted a large number of imported cases from Dhaka. All the districts that reported five or more probable cases (Chittagong, Bogra, Dhaka, Munshiganj, Narshingdi) to IEDCR were amongst our top twenty districts based on the predicted cumulative number of importations from Dhaka. Our model also predicted importations to areas far from Dhaka, such as Joypurhat district, whereas the naive diffusion model predicted no importations. Cases from districts outside Dhaka were only reported for a short period of time, between July 17 and August 10, and thus, the start of the outbreak could not be inferred from the data to compare with our importation estimates. Nonetheless, the timing of cases is consistent with increased likelihood of importation during Eid, and the long serial interval (≈ 23 days [46]) for chikungunya.