Malaria data
The malaria incidence data (number of recorded cases per 1000 persons per year) were collected by the National Malaria Elimination Programme from January 2015 to August 20187 and the population sizes were from the census in 201113.
Human movements
We estimated human movements between all pairs of unions (smallest administrative unit in Bangladesh) with at least one mobile phone tower using mobile phone call data records from April 1st to September 30th, 20177. Subscribers were assigned locations based on the union where their most frequently contacted mobile phone tower was located on a given day. Movements were identified if the primary location changed from 1 day to the next. We first calculated the average numbers of daily movements among unions, then obtained the proportion of trips from union i to union j, denoted as pij, by dividing the average number of daily trips from union i to j by the total number of daily trips originating from union i. Since human movements to and from unions without any mobile phone towers could not be estimated from the mobile phone calling data—these unions are primarily located in forested areas with higher malaria incidence rates—we trained a statistical model based on estimated human movements from the mobile calling data to estimate human movements for unions without any mobile phone towers. We used a modified version of the gravity model, fit to observed data for locations with mobile phone towers, to estimate the movement between unions that were missing towers, Mij. Specifically, we assumed Mij ~ Poisson(λij) and fit the following Poisson regression model:
$$log ({\lambda }_{{{{\rm{ij}}}}})={\alpha }_{{{{\rm{i}}}}}+{\alpha }_{{{{\rm{j}}}}}+{\beta }_{1}{Log}({H}_{{{{\rm{i}}}}})+{\beta }_{2}{Log}({H}_{{{{\rm{j}}}}})+{\beta }_{3}{T_{ij}}$$
(1)
where λij is the expected number of trips from union i to union j; αi and αj are fixed-effects at the district-level for unions i and j, respectively; Hi and Hj are the population sizes of unions i and union j, respectively; and Tij is the travel time between union i and union j14. We used the glm package in R to estimate the parameters, and used the fitted model to estimate the unobserved Mij.
We also estimated the number of subscribers who remain in the same union, i.e. the number who do not move in a given day, based on a similar model fit to observed data for locations with mobile phone towers. We added income as a covariate to the model and fit the following Poisson regression model:
$$log ({\lambda }_{{{{\rm{ii}}}}})=\alpha +{\beta }_{1}{Log}({H}_{{{{\rm{i}}}}})+{\beta }_{2}{{Inc}}_{{{{\rm{i}}}}}$$
(2)
where λii is the expected number of subscribers who remain in a given union; α is the fixed-effect at the district-level for union i; Hi is the population size of union i, as defined in Eq. 1; and Inci is the average income of union i from15. Missing pij values were then calculated according to \(\frac{{\lambda }_{{ij}}}{{\sum }_{j}{\lambda }_{{ij}}}.\)
Metapopulation malaria transmission model
We used a metapopulation malaria transmission model, developed by Ruktanonchai et al.8, parameterized by the human mobility data described above. This model is an extension of the Ross-Macdonald framework, which has been adapted to include human movement dynamics8,16,17. Within this model, we track the proportions of infected humans (Xi) and mosquitoes (Yi) for each union i.
The model incorporates two main aspects of human mobility that affect malaria transmission. Firstly, it considers the impact of infected human travelers on the infection of susceptible mosquitoes in their destination. This is quantified by the variable \({\kappa }_{i}\), which represents the proportion of infected humans in location i, including both residents and visitors. Secondly, the model accounts for how infections of residents in each location i are influenced by exposure to infectious mosquitoes both locally and in locations they visit. This influence is quantified by aggregating the impacts from all locations, with each contribution weighted by the mobility estimates pij.
$${\kappa }_{i}=\frac{{\sum }_{j}{p}_{{ji}}{X}_{j}{H}_{j}}{{\sum }_{j}{p}_{{ji}}{H}_{j}}$$
(3)
$$\frac{d{X}_{i}}{{dt}}={\sum }_{j}{p}_{{ij}}{m}_{j}{ab}{Y}_{j}\left(1-{X}_{i}\right)-{{rX}}_{i}$$
(4)
$$\frac{d{Y}_{i}}{{dt}}={ac}{\kappa }_{i}\left({e}^{-\mu \tau }-{Y}_{i}\right)-\mu {Y}_{i}$$
(5)
Here, m represents the ratio of the total female mosquito population to the total human population, r describes the rate at which infected humans recover, μ indicates the mortality rate of infected mosquitoes, and τ refers to the incubation period of the disease within mosquitoes. The biting rate of mosquitoes on humans is denoted by a, whereas b and c represent the probabilities that a bite from an infectious mosquito will successfully transmit the disease to a susceptible human and vice versa, respectively. Hi denotes the human population size in union i13. Among these parameters, m is more likely to be study-site specific, while other parameters representing features of the mosquito biology and malaria infection and transmission and are commonly assumed to be constant in previous studies8. Therefore, we used parameter values for a, b, c, r, μ, and τ from other malaria studies (see Supplementary Table 1) and solved for m for each union using the following approach.
Since mosquito dynamics are relatively faster compared to human dynamics, we assumed a quasi-equilibrium for infectious mosquitoes. We solved for quasi-equilibrium Yi by setting Eq. 5 to zero, and substituted this value into Eq. 4 to derive the resulting equation for the proportion of infectious humans (Eq. 6) (see details in ref. 8). The solution for quasi-equilibrium Yi, Yi*, is \(\frac{{ac}{\kappa }_{i}}{\mu +{ac}{\kappa }_{i}}{e}^{-\mu \tau }\).
$$\frac{d{X}_{i}}{{dt}}={\sum }_{j}{p}_{{ij}}{m}_{j}{ab}{Y}_{j}^{* }\left(1-{X}_{i}\right)-{{rX}}_{i}$$
(6)
Equation 6 describes how the proportions of infected humans change over time. We followed the methods developed by Ruktanonchai et al.8, and described briefly below, to estimate the vectoral capacity, mi, for each location. To solve for mi, we assume steady-state and set Eq. 6 to zero.
\({\sum }_{j}{p}_{{ij}}{m}_{j}{ab}{Y}_{j}^{* }\left(1-{X}_{i}\right)-{{rX}}_{i}=0{\to }^{{\mbox{yields}}}{\sum }_{j}{p}_{{ij}}{m}_{j}{ab}{Y}_{j}^{* }=\frac{{{rX}}_{i}}{1-{X}_{i}}\), which can be expressed in matrix form as
where
$$A=P{{{\rm{diag}}}}(\, f\left(X\right))$$
$$O=\left(\begin{array}{c}{O}_{1}\\ \vdots \\ {O}_{n}\end{array}\right)\; {{{\rm{with}}}}\; {O}_{i}=\frac{{m}_{i}{a}^{2}{e}^{-\mu \tau }}{\tau },$$
$$g\left(X\right)=\left(\begin{array}{c}{g}_{1}\left(X\right)\\ \vdots \\ {g}_{n}\left(X\right)\end{array}\right)\; {{{\rm{with}}}}\; {g}_{i}\left(X\right)=\frac{r{X}_{i}}{1-{X}_{i}}\; {{{\rm{and}}}}$$
$$f\left(X\right)=\left(\begin{array}{c}{f}_{1}\left(X\right)\\ \vdots \\ {f}_{n}\left(X\right)\end{array}\right) \; {{{\rm{with}}}}\; {f}_{i}\left(X\right)=\frac{{bc}{\kappa }_{i}\mu }{{ac}{\kappa }_{i}+\mu }.$$
Then O can be solved by
and
$${m}_{i}=\frac{{O}_{i}\tau }{{a}^{2}{e}^{-\mu \tau }}.$$
Since our observed data is malaria incidence by location, Ii, we used the steady-state relationship between incidence and prevalence and set \({X}_{i}^{* }=\frac{{I}_{i}}{r}\). Thus, given the malaria incidence of each union, we calculated analytical solutions for mi by solving for O (Eq. 8).
This analytical steady-state solution does not restrict values of mi to be >zero. For regions with very low or zero incidence, it is possible for the solution of mi to be negative. However, as negative mi has no biological meaning, for the unions with derived mi −10). To ensure this did not influence the fitting to the incidence data, we compared the incidence predicted from this mechanistic model and the empirical incidence data, and found they were consistent (Supplementary Fig. 1).
For comparison, we also constructed the basic Ross-Macdonald model16 without spatial component as follows:
$$\frac{d{X}_{i}}{{dt}}={m}_{i}{ab}{Y}_{i}\left(1-{X}_{i}\right)-{{rX}}_{i}$$
(9)
$$\frac{d{Y}_{i}}{{dt}}={ac}{X}_{i}\left({e}^{-\mu \tau }-{Y}_{i}\right)-\mu {Y}_{i}$$
(10)
Similarly, we assumed quasi-equilibrium for Yi by setting Eq. 10 to zero, substituting this value into Eq. 9, and then setting the modified Eq. 9 to zero to derive the expression for mnomob,i without considering spatial dynamics:
$${m}_{{{{\rm{nomob}}}},\,i}=\frac{r{X}_{i}({ac}{X}_{i}+\mu )}{{a}^{2}{bc}{e}^{-\mu \tau }{X}_{i}(1-{X}_{i})}$$
(11)
With two versions of mi, one considering mobility and one without, we derived two corresponding versions of R0 for each union i by
$${R}_{0i}=\frac{{m}_{i}{a}^{2}{bc}{e}^{-\mu \tau }}{r\mu }.$$
(12)
In Eq. 6, the changes in the proportion of infectious humans in location i are driven by infectious mosquitoes either in location i (\({p}_{{ii}}{m}_{i}{ab}{Y}_{j}^{* }\)) or in location j (\({p}_{{ij}}{m}_{j}{ab}{Y}_{j}^{* }\)). To compute the proportion of infected individuals whose residential location is union i and who were infected either while staying in union i (Cii) or while traveling to union j (Cij), each term is divided by the sum of all terms as follows:
$${C}_{{ij}}=\frac{{p}_{{ij}}{m}_{j}{ab}{Y}_{j}^{* }}{{\sum }_{k=1}^{n}{P}_{{ik}}{m}_{k}{ab}{Y}_{k}^{* }}$$
(13)
The number of unions is denoted by n. The proportion of imported infections (can be viewed as “sink score”) in union i is then calculated by \({\sum }_{j=1\,(j\ne i)}^{n}{C}_{{ij}}\). The “source score” for union i is equal to \({\sum }_{j=1\,(j\ne i)}^{n}{H}_{j}{I}_{j}{C}_{{ji}}\), where Ij is the incidence rate in location j.
Finally, we defined highly populated areas (Hhigh) as those with population sizes in the top quartile. We quantified the connectivity to these areas for each union by summing the number of people traveling to or from the highly populated areas as follows:
$${D}_{i}={H}_{i}{\sum }_{j\in {H}_{{{\mbox{high}}}}}{P}_{{ij}}+{\sum }_{j\in {H}_{{{\mbox{high}}}}}^{n}{H}_{j}{P}_{{ji}}$$
(14)
To identify malaria elimination strategies, we simulated the impact of local elimination on overall reduction of incidence in the region. We did so by setting the mosquito-to-human ratio for each union to zero one at a time, which is analogous to the maximum level of mosquito control. The effect of interventions was calculated by the reduction in the number of infected individuals across the whole CHT region.
Ethical approval
Ethical approval was obtained from the Oxford Tropical Research Ethical Committee (1-15), Bangladesh Medical Research Council Ethical Committee (BMRC/NREC/2013-2016/1154) and Harvard University Human Research Protection Program (IRB14-2669). No consent was required as we used anonymized routinely collected malaria surveillance data which was aggregated as numbers of cases with no personally identifiable information. Ethical approval to use this dataset was approved by all three committees.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
