Study site overview
In 2012–16, health facility catchment areas (HFCAs) in Gwembe District near Lake Kariba in Southern Province, Zambia, took part in operational-scale mass-test-and-treat (MTAT) [23] and mass drug administration (MDA) studies [24]. During these operations, data were collected on household locations, individual infection status as measured by rapid diagnostic test (RDT), recent symptoms and health-seeking behaviour, household ownership and individual usage of insecticide-treated nets (ITNs), and whether the household had recently received indoor residual spraying (IRS).
Household model
All simulations were run with EMOD v2.10 [25], a stochastic agent-based model of malaria transmission that includes vector life cycle dynamics and within-host immunity effects [26, 27]. Relationships between force of infection, incidence, prevalence, and infectiousness have been calibrated to field data from study sites in sub-Saharan Africa [28]. Each household was modelled as a separate locus of transmission; within each household, transmission was homogeneous with biting propensities modulated by human size and individual ITN usage.
Household locations and sizes (Additional file 1) were constructed from a master list of households assembled from GPS locations and household rosters which included ages, genders, and household relationship structures recorded in the 2012–13 MTAT rounds data [23].
Human movement
Each household had a randomly selected set of five other households within the HFCA that residents were permitted to visit, with each stay lasting an average of 3 days. For within-HFCA movement, each individual made on average 0.075 visits per year to other households. Migration rate within HFCAs was constant throughout the year.
In addition to moving between households within the same HFCA, individuals could also travel to an external high-transmission area that approximated lakeside regions. These lakeside trips lasted an average of 30 days. The external high-transmission area had a population of 1000 people, and residents of this area were not permitted to travel to the HFCAs. Because Munyumbwe HFCA is directly connected to lakeside areas by a major road, migration between Munyumbwe and the external high-transmission area was modelled at ten times the rate as for Bbondo and Luumbo HFCAs. In Bbondo, 95% of simulation-years had 38–65 individual trips to the external high-transmission area, with mean 0.0065 trips per person per year and 0.5–0.7% of Bbondo’s population visiting the external high-transmission area in any given year; in Luumbo, 95% of simulation-years had 67–97 trips, mean 0.0068 trips per person per year, and 0.5–0.8% of Luumbo’s population visiting the external high-transmission area in any given year; and in Munyumbwe, 95% of simulation years had 1569–1756 trips, mean 0.067 trips per person per year, and 6.3–6.9% of Munyumbwe’s population visiting the external high-transmission area in any given year. Following patterns seen in the MTAT rounds data in Munyumbwe HFCA, migration to the external high-transmission area was 10% higher between June and August. In Bbondo and Luumbo HFCAs, insufficient data was available to estimate relative seasonal migration rates, and migration rate to and from the external high-transmission area was modelled as constant throughout the year.
Vector species and movement
Based on entomological studies conducted in the Lake Kariba region and other entomological data from Zambia ([29,30,31]; personal communication from Javan Chanda), both Anopheles arabiensis and Anopheles funestus vectors were included in the models. Anopheles funestus was parameterized with 95% indoor biting and Anopheles arabiensis with 50% indoor biting to account for outdoor resting behaviour, which limits the impact of ITN and IRS effect sizes [31]. These indoor biting fractions influence vector susceptibility to vector control interventions such as ITNs and IRS. A. arabiensis was modelled with rainfall-driven habitat peaking in February and a small amount of year-round available habitat, while A. funestus was modelled with a lagged rainfall-driven habitat component peaking in April and a dry season component peaking in November to reflect marshy areas that become fast-moving streams during the rainy season.
Vector numbers were determined by tuning the abundance of larval habitat associated with each household and modulated by seasonal patterns of rainfall; this process is described below in “Transmission intensity”. Each vector was modelled as an individual agent. Existing data on anopheline flight patterns is still sparse and highly dependent on experimental conditions and local geography [32], so a simple model of vector movement was implemented. Vectors migrated between households within an HFCA with distance-dependence as in the profile shown in Additional file 2, which was chosen based on literature review and calibration to mark-release-recapture studies [33, 34]. Vectors moved between households within 500 m apart at a high but decreasing rate depending on inter-household distance and could make trips of up to 3 km with a small probability; 90% of vector flights were within 370 m and 97% within 500 m. The distance-dependence of vector migration between pairs of households was weakly modulated by a propensity to migrate toward households with abundant associated larval habitat. Vectors did not migrate between HFCA households and the external high-transmission area. While local vector and human movement patterns were not explicitly varied in this study, the impact of spatial extent of transmission was measured indirectly by testing intervention mixes in areas of different household density.
Intervention coverage
The timing and spatial coverage of ITN usage, IRS coverage, mass drug campaign coverage, and treatment-seeking rates were modelled according to self-reported data collected during mass treatment rounds (Additional files 1, 3, 4, 5, 6, 7).
ITN parameterization
In simulation, ownership of an ITN implies nightly usage; in that way, ITN “coverage” and “usage” are interchangeable terms when referring to simulations. The fraction of people reporting having slept under a net last night during the drug campaign rounds was used to parameterize ITN coverage rates. To determine local ITN coverage, each HFCA was gridded at 1 km, ITN coverage was calculated for each grid cell, and household ITN coverage was modelled at the level seen over the grid cell (Additional files 3, 4, 5). As observed during campaign rounds, young children and adults were given a higher likelihood of receiving a net than school-age children [21]. At each MTAT or MDA/fMDA round, individuals were identified who reported sleeping under a net they received after the previous round of data collection. These nets were then distributed according to the reported net age and gridded coverage. ITN efficacy was parameterized with initial efficacy of blocking and killing of 0.9 and 0.6 and half-lives of 2.5 and 1.5 years, respectively for both vector species [35]. Based on net usage and ages reported across campaign rounds, individuals were modelled to discard ITNs after a mean retention time of 9 months. ITN distributions after 2016 were modelled as occurring in on June 15.
IRS parameterization
Indoor residual spraying coverage was parameterized with a similar approach to ITNs. IRS timing was backdated from rounds data reporting on how long ago houses were sprayed. Campaigns after the Jan 2016 spraying were assumed to continue as planned and were modelled occurring in early January. Each household within a grid cell received IRS with probability equal to the fraction of households in that grid cell reporting having received IRS that year (Additional file 6). IRS was parameterized similarly for both vector species. Prior to 2015, IRS had initial killing efficacy of 0.5 and half-life of 2 months. Beginning with the Jan 2015 spray campaign, spraying was modelled with initial killing efficacy of 0.6 and half-life of 18 months to reflect the switch to Actellic, a longer-lasting insecticide; this parameterization results in 20–25% reduction in vector numbers in the subsequent year when 50% of households receive IRS. Since Bbondo HFCA consistently had little to no reported spraying, with at most 4.5% of households reporting recent spraying in any round, IRS campaigns were not modelled there.
Mass drug campaigns
Mass test-and-treat was distributed in June, August, and October of each of 2012 and 2013, and a pilot round of MTAT in December 2011 was also included for Bbondo and Munyumbwe HFCAs. Drug campaign coverage achieved in each HFCA during the MTAT rounds was estimated from the fraction of longitudinally linked individuals found within that HFCA in subsequent campaign rounds, corrected by the estimated linkage efficiency [21]. During simulation, MTAT and MDA coverage is random at the individual level, not correlated between rounds or within households. In Bbondo and Munyumbwe HFCAs, MTAT coverage was modelled as 60%, except for the pilot round, which was set to 30 and 40% respectively for Bbondo and Munyumbwe. In Luumbo HFCA, MTAT coverage was modelled at 40%, except for the June 2012 round, which was modelled at 20%. During MTAT, all covered individuals who tested positive by RDT received a full course of artemether–lumefantrine (AL). RDT sensitivity was assumed to be 40 parasites per µL when administered by campaign staff.
For the 2014–2016 mass drug campaigns, coverage was assumed to be the same as the MTAT coverage: 60% for Bbondo, 40% for Luumbo, and 60% for Munyumbwe. Bbondo and Luumbo HFCAs received MDA, where all individuals received a full course of DHA–piperaquine (DP) regardless of infection status, and Munyumbwe HFCA received focal MDA (fMDA), where only individuals residing in a household with someone who tested positive by RDT received DP [24]. MDA or fMDA rounds were distributed in December 2014, February 2015, September 2015, and February 2016.
Case management
Simulated case management rates prior to 2014 reflected campaign data collected during the 2012−2013 MTAT rounds, which showed both age- and distance-dependence of health-seeking behaviour [21]. Children under 15 years of age had a relative 50% greater chance of receiving care than adults. All clinical and severe cases that received treatment received an age-appropriate dose of AL within 3 days of presenting with symptoms.
Baseline case management rate modelled health-seeking behaviour depending on distance from the local clinic as observed during the MTAT rounds (Additional file 7). Case management rates after the CHW programme began in mid-2014 were assumed to increase as follows: cases in households within 1 km of the CHW’s location sought treatment at 60% for adult cases and 90% for child cases, and households beyond 1 km of a CHW followed the same distance-dependence as was observed for the clinics during the MTAT rounds. CHW and clinic locations were obtained from Zambia’s DHIS2 database [36].
Reactive case detection
During RCD, programme guidelines indicate that a treated clinical case should trigger a CHW to test-and-treat all individuals in households within 140 m of the index case [8]. This was replicated in the model. Sensitivity of RDTs given during RCD was assumed to be 100 parasites per µL, a lower sensitivity than during MTAT as CHWs may have more difficulty with proper interpretation of RDT results than campaign staff [37], and all treatments used AL.
Since CHWs did not report which active follow-ups are triggered by which passive cases, it is difficult to ascertain the rate at which index cases were followed up with reactive activities and whether the follow-up procedure correctly followed the 140 m-radius guidelines. During follow-up activities, standard RCD was parameterized as 40% of index cases triggered follow-up, test-and-treat conducted within 140 m of index household, and 60% of people available for test-and-treat during follow-up [11]. Solely for out-of-sample comparisons in Munyumbwe HFCA, RCD was modelled as follows: in the northeast area where there is no CHW, no RCD was modelled. In the southwest region where transmission has historically been very high, 20% of index cases triggered follow-up, follow-up occurred only within the household of the index case, and 40% of people were available for test-and-treat. In the central region where the major roads intersect and many CHWs are present, 50% of index cases triggered follow-up, follow-up occurred within 50 m of the index household, and 40% of people were available for test-and-treat.
In hypothetical intervention scenarios, RCD was always modelled with a radius of 140 m from the index case. “Perfect” RCD indicates 100% of index cases occurring in the HFCA triggered follow-up activities and 100% of people were available for test-and-treat during follow-up. In all RCD scenarios, follow-up activities were modelled to occur the same day as treatment of the index case. RCD distributed with focal test-and-treat (fTAT) gave drugs only to individuals testing positive by RDT while RCD distributed with focal drug administration (fDA) gave drugs to all individuals found within the target radius. In all simulations, no upper limit on the number of people tested or treated with RCD was imposed.
Transmission intensity
Household-associated larval habitat abundances and intervention coverages together determine each household’s baseline transmission intensity. Each household’s habitat abundance contributed to local vector numbers, household biting rate, and thereby individual risk of testing positive by RDT. Larval habitat abundances were selected from 200 sets of random household habitat abundances by comparing resulting simulation outputs to two types of data from the 2012–13 MTAT rounds: (1) HFCA-wide prevalence by RDT measured during the seven rounds of MTAT between 2011 and 2013, and (2) from the June 2012 data, the probability of an individual testing positive given another positive individual within the same household, within 50 m but not in the same household, and between 50 and 200 m away (Additional file 8). In addition to HFCA-wide prevalence by RDT, prevalence measurements for subregions in Luumbo and Munyumbwe were also considered as calibration targets as village-scale baseline prevalence in Luumbo and Munyumbwe showed considerable heterogeneity (Additional file 9).
Comparisons between simulations and observed data were calculated using euclidean distance, and distances for all comparisons were summed to calculate overall score. The 20 lowest-score habitat abundance samples were used for out-of-sample predictions and scenario projections. Data collected during the 2014–16 MDA and fMDA campaigns, routine clinical case counts from health facilities, and data collected by CHWs were not used for model calibration.
In Bbondo HFCA, all households had the same ratio of arabiensis and funestus habitat availability. In Luumbo and Munyumbwe, households within each subregion had the same ratio of arabiensis and funestus habitats; subregion-specific ratios were chosen based on preliminary entomological data in nearby sites to reflect subregion-specific seasonality.
Scenario projections
For all HFCA models, all vector control (ITN and IRS) and mass drug distributions (MTAT, fMDA, MDA) interventions prior to 2017 were included in the intervention scenario packages (Fig. 2). Various potential case management rates, RCD coverages and methodologies, vector control distributions, additional mass drug campaigns, and malaria control measures implemented in the external high-transmission area were simulated over the period January 2014 through December 2021.
All intervention packages in the high-transmission area included annual IRS at 50% coverage and 30% case management rate. In addition, each intervention package included: (package 1) ITNs distributed every 3 years at 30% coverage; (package 2) ITNs every 3 years at 60% coverage; (package 3) ITNs every 3 years at 80% coverage; (package 4) ITNs every 2 years at 80% coverage; (package 5) ITNs every 2 years at 80% coverage and two MDA rounds each year at 70% coverage.
Each intervention scenario was simulated with 1000 stochastic realizations, 50 for each of the 20 best habitat abundance samples, to capture the range of possible outcomes. Realizations of each scenario where zero locally-acquired infections occurred over the entire 2020 calendar year were counted to have achieved elimination [38]. Locally-acquired infections included all new infections acquired in Bbondo, Luumbo, or Munyumbwe HFCA, including completely asymptomatic infections, but not infections acquired in the external high-transmission area and imported into Bbondo, Luumbo, or Munyumbwe.