Data sources
Much of the data used to design and calibrate the simulation model of malaria transmission on Bioko Island was collected by the BIMEP and the National Malaria Control Programme (NMCP). Population census data were collected as part of two campaigns for distributing long-lasting insecticidal nets across the island in 2015 and 2018. The population present during these visits was counted. On average, approximately 12% of the houses had no response or their households denied access to distributors and could not be included in the census (Fig. 1b) [20,21,22]. Hence, these census data included approximately 88% of all households on the island and, despite the marginal under-count, are regarded as the most accurate estimate for Bioko and used as the denominator for all interventions undertaken by the BIMEP [23]. (To clarify, bed nets are also made available to all households in all years on Bioko Island, not only to the households counted in the surveys). Since the beginning of malaria control on Bioko Island in 2004, the BIMEP has performed extensive annual malaria indicator surveys (MIS) that have collected epidemiological, demographic, and socioeconomic data. The present study draws upon MIS data collected in each of the years from 2015 to 2018, which sampled an average of 16,500 respondents each year, to analyse transmission patterns on Bioko Island [24,25,26,27]. MIS data were again collected in 2019 and 2020, but the analyses in the present manuscript were performed before those data were available.
The MIS data include results from rapid diagnostic tests (RDT; CareStart Malaria Pf/PAN (HRP2/pLDH) Ag Combo RDT, AccessBio Inc, Monmouth, USA), used to diagnose survey participants for the presence of detectable Plasmodium falciparum parasites (Fig. 1c). The RDT results in the MIS data make it possible to map the spatial distribution of occurrences of malaria infection. Guerra et al. utilized the RDT data to produce geostatistical estimates of P. falciparum parasite rate (PfPR) [17]. The geostatistical estimation techniques use the RDT results together with geospatial environmental covariates to estimate malaria prevalence across the island, accounting for the uneven distribution of samples taken by the survey [28]. These estimates reveal how the prevalence varies across geographical space, with the highest prevalence on the island occurring along the northwest coast and in the southeast; a reduced prevalence in the high population density, urbanized areas in the capital city of Malabo; and reduced prevalence in the elevated regions nearer to the centre of the island. The overall mean estimated PfPR is around 0.12 for the entire island. Note that the PfPR estimates reflect the current level of interventions on the island. The PfPR estimates are made for all age groups, as there were no significant differences between PfPR measurements made for children (ages 2–10 years) and the overall sample population [17]. The median surface of PfPR estimates are shown in Fig. 1d, reproduced with permission from [17].
The MIS asks respondents about recent travel history, providing information on where it was that people could have become exposed. For the MIS in 2015–2017, respondents reported whether they had recently taken a trip where they had spent at least one night away from their home residence in the preceding 8 weeks. Respondents who reported at least one trip also reported one of seven possible destinations, to each of the regions of Malabo, Baney, Luba, Riaba, Moka, Ureka, or off-island (Fig. 1a). Most of the travellers (84%) reporting off-island travel reported going to mainland Equatorial Guinea [17]. The 2018 MIS included an additional question on travel duration, making it possible to assess how many days travellers had been away from home [29].
Given that so many travelers reported going to mainland Equatorial Guinea, the model also required knowing estimates of prevalence in that region. The Malaria Atlas Project estimated the overall median PfPR in the region of mainland Equatorial Guinea to be 0.43 in children aged 2–10 at the time of MIS data collection [15]. This estimate is corroborated by a study conducted in 2015 by Ncogo et al., which reported an overall PfPR in the region to be 0.46 [16] among all age groups. The study by Ncogo et al. did find geographical variation in PfPR, with elevated prevalence as high as 0.58 in rural areas and 0.34 in urban areas [16], but the MIS travel data do not contain enough detail to report exactly where travellers spent their time in the mainland. For this reason, along with the fact that even in lower-prevalence areas the PfPR remains far higher than on Bioko Island, the median PfPR estimate is used for the mainland region. That the prevalence is so much higher on the mainland than on Bioko Island suggests that travellers would experience a much higher risk of infection while visiting the mainland than they would at home. This is consistent with the prior studies and analyses performed which suggest that there is an elevated risk of prevalence among those who had recently travelled to mainland Equatorial Guinea [6, 9, 17].
Lastly, the MIS data also included some basic information on treatment-seeking behaviour. Respondents reported whether they had recently experienced a fever and whether they sought treatment for malaria.
Designing the simulation model
A family of models was constructed in order to simulate malaria transmission patterns on Bioko Island, based on the census and epidemiological data collected from 2015 to 2018. The models are continuous-time, event-driven agent-based stochastic simulations that describe how individual human hosts become infected through contact with mosquito vectors in the environment; how malaria infections run their course; and how infected individuals contribute to onward transmission of subsequent infections. The models are spatially explicit, meaning that prevalence and local transmission intensity are allowed to vary across geographical space. The model also simulates human travel behaviour, and human hosts may experience different levels of transmission risk as they move from one location to another. Parameters describing local transmission and human mobility were fit to be consistent with the geospatial analysis of prevalence and MIS data (see below). Code supporting the simulations may be found in the macro.pfsi directory at https://github.com/dtcitron/bioko_island_travel_materials, and an explanation for how the simulation program works may be found in Additional file 1: “Simulation software description” section.
Within the simulation, the resident population of Bioko Island is based on the 2018 census data, which were aggregated into 241 gridded 1 km × 1 km map-areas [22] (as shown in Fig. 1b). Each populated map-area from the census is represented as an isolated patch in the simulation, home to as many human hosts as in the corresponding map-area from the census data. Patch residents spend most of their time in their home patches but occasionally travel to other ones. To simulate malaria importation, a 242nd patch is included to represent the destination of off-island travel (mainland Equatorial Guinea). This final patch serves as a boundary condition for the simulation: Bioko Island residents import cases through experiencing exposure risk while traveling off-island, and the 242nd patch represents the off-island transmission environment.
The core of the simulation model represents transmission of malaria parasites between human hosts and the vectors. The core of the transmission model is based on the Ross–Macdonald model, which includes a modeled description of human hosts, the population dynamics of the vector population, how the vectors and humans interact with one another, and the course of infection in a malaria-afflicted individual [30, 31].
Figure 2 shows a simplified schematic representing the modelled compartments for both mosquitoes and humans within a single patch. The model tracks mosquito population dynamics: within each patch, adult mosquitoes emerge at a rate which reflects the local ecology and die at a constant rate. Adult mosquitoes become infected with malaria parasites when they blood feed on infectious human hosts. After the extrinsic incubation period, surviving adult mosquitoes then become infectious. Each day, within each patch a number of infectious bites is generated based on local population of infectious mosquitoes and distribute those infectious bites across the human hosts who are present that day. Thus, susceptible human hosts can become infected. Infected human hosts may develop symptoms and consequently seek treatment, at which point their infections become cleared and they remain protected against new infections for a short period. Two final model features are not shown in Fig. 2. The complete transmission model simulates together many patches, and allows for human hosts to travel between them. The simulation model also has the capability to simulate the distribution of a pre-erythrocytic vaccine, where vaccinated people experience a lower infection risk.
Calibrating the simulation model
The key to calibrating each simulation model is to set the mosquito population density such that it sustains a level of transmission intensity, which subsequently produces the correct PfPR within each patch. The simulation is calibrated using a geospatial estimate of PfPR within each patch drawn from the Bayesian posterior probability distribution [15, 17], and then the transmission model is used to translate PfPR into a local force of infection (FOI).
Each patch’s FOI is a quantitative representation of the transmission risk experienced by human hosts found there. Accounting for human movement is a key part of this analysis, and the procedure for calculating FOI in each patch is adapted from the source-sink analysis described in [32, 33] (refer to Additional file 1: “Model calibration” section for more details). From the perspective of an individual human host, their risk of becoming infected with malaria is the result of the FOI they experience. The average FOI is computed as a sum of the FOI experienced in each patch visited by that individual weighted by the duration of time spent in each of those patches.
Knowing the FOI, the transmission model is again used to derive that patch’s mosquito population density. Within the simulation, the mosquito population density produces blood feeding activity, which in turn affects the transmission risk (FOI) which finally gives rise to the geospatial PfPR estimates when the mosquito and human populations are coupled together within the simulation. Thus, within the simulation it is possible to calibrate the vector ecology across the different patches on the island such that it reproduces the geospatial PfPR estimates.
Incorporating travel data
Accounting for human hosts’ travel patterns is the last step required for properly calibrating FOI. The total FOI experienced by an individual host includes both FOI at home as well as FOI experienced while travelling away from home. For this, the MIS travel data are used to construct and parameterize a model of human travel patterns on Bioko Island [29]. Each human host’s travel behaviour is modeled in three steps: the human host chooses when to leave their home patch; chooses the destination; and chooses how many days they spend away before returning to their home patch. Within the simulation, each individual requires a set of parameters that includes the frequency of travel; the multinomial probability distribution for determining travel destination; and the mean duration of their trip [34]. This model of movement behaviour is too simple to allow for trips with multiple destinations, but it is nevertheless consistent with the MIS travel history data. The frequency of leaving home is estimated based on the frequency of trips reported by MIS respondents. The probability of travelling to each travel destination is estimated based on the relative frequencies of trips from each patch to each destination region. The duration of trips is estimated based on the distribution of trip duration reported—an average of about 10.5 days for travel within Bioko Island and 20 days for travel to mainland Equatorial Guinea. Data for fitting travel frequencies and probabilities were available from the MIS from 2015 to 2018, but data for fitting travel duration were only available with the most recent 2018 MIS [29]. Refer to Additional file 1: “Movement model parameterization from data” section for a detailed description of the movement model’s parameterization. For simplicity, all human hosts from each patch have the same set of movement parameters, where those parameters are allowed to vary from patch to patch according to the MIS data. For example, individuals who live in the southern parts of Bioko Island travel more frequently and tend to choose Malabo as the destination, whereas individuals who live in Malabo travel less frequently but tend to choose mainland Equatorial Guinea as the destination. The model ignores migration of mosquitoes between patches.
Estimating uncertainty
One of the benefits of using a simulation model is that it enables producing results which reflect the uncertainties underlying the data used to calibrate the model. The geostatistical PfPR estimates include mean surfaces across the island (plotted in Fig. 1d) as well as a full ensemble of draws from the joint posterior distribution. Each draw from the joint posterior distribution by itself represents a map of PfPR estimates across Bioko Island that also accounts for spatial correlations between the PfPR in different patches. The full ensemble of draws from the joint posterior distribution represents a sample which reflects the uncertainty of the mapped PfPR estimates [15]. An ensemble of simulated outcomes is constructed using the ensemble of draws. For each simulated scenario generates1000 simulation runs. Each simulation run is calibrated using its own PfPR surface draw. The full ensemble of simulation runs represents a family of models whose outputs reflect the variability in simulated outcomes when considered together. These variations are expressed in terms of error bars or confidence intervals. The error bars shown in the subsequent plots reflect both the uncertainty inherent to the stochastic simulation as well as the uncertainty associated with the PfPR maps.