Overview
A three-step process was used to generate LGA-level predictions of potential national strategic plans (Fig. 1). At the outset, the goal was to capture the intrinsic potential of each LGA to support malaria transmission in a baseline period before 2010 when most interventions were not scaled up nationwide. Data and geospatial modeled surfaces from 2010 or before were used to group LGAs into epidemiological archetypes. For each archetype, baseline malaria transmission was calibrated to 2010 data. Next, Nigeria’s intervention history from 2010-20 at the LGA level was imposed on the baseline models to generate 774 LGA-level models up through 2020. Last, various future intervention strategies were applied to the LGA models and intervention impact on prevalence, incidence, and mortality was assessed. The following sections summarize the study methods. Additional details are provided in the Supplement.
Clustering LGAs with similar baseline transmission patterns to generate epidemiological archetypes
Calibrating baseline transmission intensity and computing intervention coverages for each LGA in the model was challenging due to the computational power required to simulate 774 LGAs, some of which lacked subnational data. To address this problem, LGAs with similar baseline malaria transmission patterns (climatic patterns, vectors, and baseline transmission intensity) were clustered (Fig. 2). The following LGA-level features were used for clustering: average monthly rainfall [22], average monthly temperature suitability index (TSI) [23], overall relative abundance of three vector species [24], annual modeled Plasmodium falciparum parasite prevalence (PfPR) values from the Malaria Atlas Project (MAP) [25, 26] for 2000, 2004, 2006, 2008, and 2010, and estimated ITN use between 2008–2010 from MAP [27]. LGAs that shared similar climatic, vector, and transmission features were grouped into 60 epidemiological archetypes using the CLARA algorithm (see Supplement) [28, 29]. Given that data from the 2010 Malaria Indicator Survey (MIS) was used to calibrate parasite prevalence, archetypes where the 2010 MIS sample population of children under the age of five years (U5) was less than 50 were reassigned to the next closest archetype (see Supplement). This final reassignment resulted in 22 epidemiological archetypes with between 13 to 92 LGAs per archetype. See Supplement for the archetype assignment of each LGA.
EMOD malaria transmission model
Malaria transmission and intervention impact was simulated using EMOD v2.20, an agent-based model of Pf transmission that comprises a model of temperature-dependent vector lifecycle and vector population dynamics, coupled to a model of human disease and immunity, and intervention effects [21, 30–33]. Parasite dynamics, immune acquisition, and clinical incidence by age were previously calibrated to field data from nine study sites in sub-Saharan Africa [34]. LGA daily air temperatures were computed using centroid longitude and latitude and data from Global Surface Summary of the Day [35]. Three local vector populations were simulated – Anopheles gambiae s.s., Anopheles arabiensis, and Anopheles funestus – with species-specific parameters for anthropophily [35–38] and probability of feeding indoors [39–41]. The number of mosquito bites per person (biting risk) was modeled with an exponential distribution for exposure [42] and surface area dependence by age to capture the heterogeneity of exposure between individuals per LGA. Population size for each LGA was rescaled to 1000 to reduce computational burden. Uniform birth and death rates were used for all LGAs [43]. Modeled LGAs were not spatially connected, and 10 infections were imported into each LGA per year.
Model inputs and scripts used for running simulations and post-processing of results are publicly available (see Data Availability section), including instructions for downloading dependencies. Scripts were written using Python 3.6.0 [44] and R 4.02 [45].
Model Calibration
i. Setting model seasonality: parameterizing monthly larval habitat availability
Seasonality of malaria within each archetype was calibrated to routine reporting of all-age uncomplicated malaria cases collected for a Rapid Impact Assessment (RIA) performed in 2019 by the NMEP. The RIA dataset contains quality-assured monthly confirmed plus presumed malaria cases from 917 nationally representative public health facilities in 601 LGAs in Nigeria for the period of 2014-18. Confirmed cases were diagnosed with either a rapid diagnostic test (RDT) or microscopy. All cases in each month were summed for each LGA and by year. The average number of monthly cases and 95% confidence intervals across all years were computed by LGA and summed by archetype. The monthly case series were rescaled such that the number of cases per month was roughly between 20 and 100 prior to attempting the simulation fit, since the model was well-behaved in this range of monthly treated cases and the extent to which the RIA data captured the true clinical burden of malaria in the archetype was unknown. For each of the 22 archetypes, seasonal variation in transmission intensity was captured in EMOD by selecting overall monthly Anopheles larval habitat availabilities such that the resulting monthly mean clinical incidence matched the shape of the monthly clinical incidence in the RIA dataset (Fig. 3A). The same seasonality was assumed for all three vector species.
ii. Accounting for non-malarial fevers
Passive surveillance, such as the RIA dataset, also contains individuals who test positive for malaria but whose symptoms are caused by a non-malarial illness, such as when an individual with an asymptomatic malaria infection becomes co-infected with another pathogen. These malaria-attributed, but not malaria-caused, fevers were accounted for in the simulations by testing a random 0.38% of the simulated population daily with an HRP2-based RDT [46], approximating a baseline rate of non-malarial febrile illness. Individuals testing positive received antimalarial treatment and were added to the daily tally of malaria cases.
iii. Setting baseline transmission intensity
Simulated baseline transmission intensity was defined as Pf prevalence detected by microscopy with a sensitivity of 50 parasites per microliter [47]. Accurately modeling parasite prevalence for each archetype in the pre-2010 model required applying a scaling factor on the monthly vector larval habitat availability to reproduce the 2010 MIS PfPR in U5, in the presence of the observed 2010 ITN usage and case management (CM) coverage (Fig. 3B and 3C). PfPR, ITN and CM coverage were computed from the MIS by taking a cluster-weighted averages across all LGAs within each archetype (clusters in the MIS and DHS survey are enumeration areas as defined by the census or an existing sampling frame). Since MIS data collection spanned a 2 — 3 month window, PfPR weighted averages were computed by archetype and month of data collection.
Fifty larval habitat scale factors were sampled for each of the 22 archetypes, resulting in 1100 unique scale factor archetype combinations. For each of these combinations, a 50-year initialization phase was run to establish population immunity in the absence of interventions. For each archetype, an ITN distribution with archetype-specific coverage as observed in MIS was applied to each of the partially immune populations generated under each of the larval habitat scale factors. The U5 PfPR during the months corresponding to the 2010 MIS sampling in the archetype was measured in each simulation, generating for each archetype a series of PfPR measurements for each larval habitat scale factor. The simulated PfPR measurements were compared to measured PfPR for the 2010 MIS using a beta binomial likelihood function and the scale factor resulting in the highest log likelihood was chosen. Each LGA was assigned the seasonality profile and larval habitat scale factor of its archetype, upon which an LGA-specific set of interventions covering the period 2010-19 was imposed.
Parameterizing historical interventions (2010–2019)
Historical interventions were parameterized for the intense intervention period from 2010–2019 when many interventions were scaled up in Nigeria (Fig. 4). Key historical intervention characteristics captured in the model were intervention type, timing, coverage, and effect size. Simulated intervention types were CM for uncomplicated and severe malaria, ITNs through household mass distribution and antenatal distributions, seasonal malaria chemoprevention (SMC), and intermittent preventive treatment in pregnancy (IPTp). Indoor residual spraying (IRS) was not included in the model since spraying activities were mostly trials or pilots covering a small subset of households [5, 48–50].
i. Case Management (CM)
The 2010, 2013, 2015, and 2018 DHS/MIS were used to estimate time-varying effective CM coverage for uncomplicated malaria by LGA (Fig. 5A). CM coverage was calculated as the fraction of children under 5 with fever in the two-week period prior to the survey treated with an artemisinin-based combination therapy (ACT). Since the DHS does not collect CM coverage data for adults, we assumed the same for adults as for children based on a 2013 study where patent medicine vendors were roughly equally as likely to treat adults and children exhibiting malaria symptoms with an ACT [51]. CM coverage estimates at the LGA level were weighted using the cluster level weights, provided within the DHS dataset. Missing CM data per LGA were replaced with the archetype-level estimate, and CM during years between DHS/MIS surveys or after 2018 used values from the most recent survey. Effective CM coverage for severe malaria was held constant at 49% for all LGAs, based on reports that injectable artesunate is insufficiently prescribed [52], not readily available in health facilities [53], and may not be sufficient in the absence of additional supportive care [54]. For almost all modeled LGAs, CM for uncomplicated malaria was lower than for severe malaria with the exception of 33 LGAs in 2018. Treated malaria cases in the model received artemether-lumefantrine with an exponential treatment delay distribution averaging 3.3 days for uncomplicated cases and 2 days for severe cases.
ii. Insecticide-Treated Nets (ITNs)
ITN coverages by age for simulated mass distributions were estimated from the 2010, 2013, 2015 and 2018 DHS/MIS surveys for the age groups 0–5, 6–9, 10–18, and > 18 years. ITN usage from the DHS/MIS is the fraction of individuals who slept under a treated bednet the night before the survey (Fig. 5B). Seasonal variation in ITN use was computed by aggregating DHS data between 2003 and 2018 and estimating monthly net use fractions. Among net users, a random 10% each night were assumed to not use their nets. ITN coverage from surveys was inflated using retention time estimates from MAP to obtain the initial ITN coverage at mass distribution, and net loss was also modeled according to half-life estimated by MAP [27]. The yearly mass distribution schedule at the LGA level was obtained from the NMEP and implemented in the simulation. However, the actual campaigns dates were assumed. LGA-level ITN coverage among pregnant women attending antenatal care in 2018 was calculated from program data shared by the NMEP and used within the simulation to distribute nets to newborns.
Nigeria’s mass campaigns have historically distributed pyrethroid ITNs. An ITN blocking rate of 0.53 was estimated from published experimental hut trials [55–57] and assumed to be uniform for these nets. ITN killing efficacy was estimated by analyzing the relationship between permethrin mortality and killing efficacy [58], which allowed parameterization by LGA and year with modeled maps of permethrin mortality for years 2005, 2010 and 2017 obtained from MAP [59] (see Supplement).
iii. Intermittent Preventive Treatment in Pregnancy (IPTp)
Dose-specific IPTp coverage in each LGA was estimated from self-reported receipt of one, two, or three or more doses of sulfadoxine-pyrimethamine (SP) in the DHS/MIS for years 2010, 2013, 2015, and 2018 (Fig. 5C). Coverage for years between surveys and up to 2020 was constructed using a monotone Hermite spline [60] (See Supplement).
Parasite prevalence outputs from the simulation were adjusted for the direct effect of IPTp among pregnant women accounting for dose dependencies. Based on Menéndez et al. [61], we estimated that one IPTp dose completely protects against infection for 10 weeks, and each subsequent dose, up to three doses, adds an additional four weeks of protection. Using these parameters, and assuming a 36-week pregnancy, we estimated the fraction of pregnancy that is unprotected for individuals who take one, two and three IPTp doses and computed adjusted prevalence estimates for pregnant women in LGA and for each simulation month. Final parasite prevalence estimates were calculated as the weighted average of the prevalence estimates in all groups, including pregnant and non-pregnant individuals (See Supplement).
iv. Seasonal Malaria Chemoprevention (SMC)
WHO recommends SMC in areas where more than 60% of transmission occurs within 4 months and annual PfPR is greater than 5% [62]. Although the LGAs receiving SMC from 2013–2017 and the number of SMC doses distributed were known, reported coverage of at least one dose exceeded 100% and the lack of a reliable population denominator made it challenging to re-estimate coverage from dose data. Post-campaign surveys were instead used to parameterize SMC coverage [63, 64]. SMC coverage was assumed to be semi-correlated between rounds, as has been observed in post-SMC surveys where the fraction of children who receive all rounds of SMC is greater than would be expected from random selection per round [8, 63]. “High access” children in the model were assumed to be more likely to receive SMC doses than “low access” children. In the simulations, children were randomly assigned to access groups at birth.
For 2013-17, 50% of U5 children were assumed to be high access and the remaining 50% low access. To capture the national SMC coverage values of between 50–80% observed in post-campaign surveys [8] (Fig. 5D), the high access group received 80–100% SMC coverage per round and the low access group received 20–60% coverage per round. Coverages in a subset of LGAs in Zamfara state was set at 100% for the high access group and 60% in the low access group based on findings from a preliminary survey [65].
For 2018-19, survey data on overall state-level SMC coverage and fraction of U5 in the high access group, defined as those who received all four SMC rounds [63, 64] for four states in 2018 and five states in 2019 was used to compute LGA-level coverage. However, some LGAs did not have low and high access children because survey data suggested that all children received the full course of SMC drugs. In those cases, children in the low access group had a coverage of 100%. The fraction in the high access group ranged from a low of 24% in Sokoto to a maximum of 53% in Zamfara. The fraction in the low access group was computed by subtracting the fraction in the high access group from the overall fraction of children (1 – fraction in high access group). LGA-level SMC coverages in the low access group did not vary within states and ranged from 20% to a high of 100%. For LGAs without survey data from 2018-19, SMC coverage from 2013-17 was used.
The pharmacokinetic parameters for the SMC drugs - sulfadoxine-pyrimethamine + amodiaquine - were obtained from the research literature [66–68] while pharmacodynamics properties were inferred through calibrations of simulated trials to clinical trial data from Ghana [69] and Tanzania [70]. Effect size of SMC on clinical incidence was validated using SMC rollout data from Mali [71]. All children were assumed to fully adhere to the 3-day treatment course.
Parameterizing effect sizes for new interventions in future scenarios (2020–2030)
i. New net formulations
The NMEP considered pyrethroid piperonyl butoxide (PBO) and Interceptor® G2 (IG2) nets for future deployment. PBO effect sizes were assumed to increase with increasing susceptibility to pyrethroids [58] and were LGA-specific, but IG2 effect was modeled as uniformly high. See Supplement for details on parameterization of kill rates and blocking rates for new nets.
ii. Intermittent Preventive Treatment in Infants (IPTi)
IPTi, historically not implemented in Nigeria, was considered in 374 eligible LGAs as part of the national strategic planning exercise for 2021–2025. Eligibility was assessed by identifying LGAs not meeting WHO recommendations for SMC and with PfPR > 10% in the absence of interventions [internal communication from WHO]. A modeled surface from MAP was used to assess PfPR eligibility. IPTi effects were applied to simulation outputs rather than incorporated dynamically. Protective efficacy parameters for IPTi were extracted from a recent systematic review [72] and multiplied by LGA-specific coverage in each scenario to generate IPTi effects on prevalence, uncomplicated cases, severe cases, and deaths. LGA-specific expected coverage for IPTi was computed by taking the mean coverage of the first, second and third doses of pentavalent diphtheria, tetanus and pertussis vaccine from the 2018 DHS survey, in anticipation of IPTi administration coupled to the expanded program on immunization for children. IPTi coverage was further increased above the expected coverage as appropriate for each intervention scenario. See Supplement for additional details.
Model validation
Simulated U5 PfPR and seasonality were validated by comparing to U5 PfPR data from the DHS and the RIA all-age routine health facility data, respectively. Prevalence comparisons were made with the Pearson correlation coefficient. State-level comparisons were made between treated uncomplicated cases, which included treated cases with non-malarial fevers in the simulation, and RIA health facility data (confirmed and suspected cases) using a cross-correlation function (CCF). The CCF at the time lag zero is a measure of the contemporaneous correlation or the linear relationship between the two timeseries. Since the RIA data was constructed from a subset of facilities of unknown catchment size, incidence values were computed using the state population as denominator and rescaled to the simulation-observed range for visual comparison. Scale factors were computed by dividing the estimated LGA monthly incidence in the simulation by the corresponding LGA incidence estimated from the simulation. The median scale factor by state and year was used to adjust monthly incidence data from the RIA. Confidence intervals were computed from a non-parametric bootstrap with 10,000 replicates using the BCa method [73, 74].
LGA intervention allocation scenarios (2020–2030)
Predictions of the impact of intervention mixes on malaria prevalence, incidence, and mortality, within each of 774 LGAs, were simulated for four scenarios (Table 1) and for 5 stochastic realizations per scenario. The NMEP directed the type, coverage, and timing of interventions to be simulated. All scenarios deployed CM, IPTp, and ITNs in every LGA, but the distribution of SMC, IPTi, and ITN type varied by scenario, as did the effective coverage of each intervention. ITN distribution schedules were every three years from the last distribution year according to NMEP data. See Availability of data and materials for link to web app with scenario maps by LGA, intervention, and year. Here, “coverage” for ITNs refers to simulated usage. The scenarios were:
-
Business-as-usual (BAU, Scenario 1). SMC distribution and coverage of all interventions (CM, standard ITNs, SMC and IPTp) remained in the same LGAs and at the same coverage levels as in 2019. All LGAs received pyrethroid ITNs.
-
National malaria strategic plan (NMSP) ramping up to 80% coverage (Scenario 2). PBO nets were distributed in all LGAs. All 404 LGAs that were eligible for SMC and 365 LGAs that were eligible for IPTi received SMC and IPTi, respectively. CM and IPTp were available in all LGAs. In Scenario 2, coverage of all interventions met or exceeded the program target of 80%. To simulate gradual intervention coverage improvement, coverage levels of CM and IPTp that were below the target of 80% in 2019 were increased over time to hit 80%, after which they were kept constant. Parameters from a beta regression model, used to estimate the average increase in CM coverage by archetype between 2013 and 2018, were scaled and added to the LGA-level 2019 coverage, until the 80% target was achieved. Coverage of ITNs, SMC, and IPTi was set to 80% for LGAs where baseline coverage (or expected coverage for IPTi) was under 80%.
-
Budget-prioritized plans submitted to the Global Fund (Scenarios 3, 4). Scenario 3 prioritized the implementation of SMC in 235 LGAs within the allocated budget, while scenario 4 included an additional 75 LGAs in the priority above allocation request (PAAR) on top of the 235 LGAs in scenario 3. Other than SMC, interventions were the same in scenarios 3 and 4. CM and IPTp coverages were assumed to increase at historical rates without exceeding 80%. Rates of increase were parameterized using the same beta regression model parameters referenced in #ii above, which was fitted to coverage data from 2013–2018 for CM and a linear regression model fitted to 2015–18 coverage data for IPTp. For LGAs where IPTp coverage had decreased in recent years, IPTp coverage was held constant for 2020-30. SMC and ITN coverages were simulated at 80%. Three net types were distributed in these scenarios: pyrethroid-only ITNs in 136 LGAs, PBO in 605 LGAs, and IG2 (α-cypermethrin and chlorfenapyr) in 33 LGAs.
Table 1
Intervention scenarios chosen by the NMEP for simulation modeling. Numbers indicate number of LGAs receiving the intervention.
| Number of LGAs targeted for each intervention |
Scenario | CM | Pyr ITN | PBO ITN | IG2 ITN | IPTp | SMC | IPTi |
1: Business-as-usual (BAU) | 774 | 774 | | | 774 | 80 | |
2: NMSP ramping up to 80% coverage | 774 | | 774 | | 774 | 404 | 365 |
3: Budget-prioritized plan with coverage increases at historical rate | 774 | 136 | 605 | 33 | 774 | 235 | |
4: Budget-prioritized plan with coverage increases at historical rate and expanded SMC in the PAAR | 774 | 136 | 605 | 33 | 774 | 310 | |
CM: Case management, Pyr ITN: Standard insecticide-treated nets/pyrethroid-only nets, PBO: Pyrethroid + piperonyl butoxide nets, IG2: (α-cypermethrin and chlorfenapyr) nets, IPTp: Intermittent treatment in pregnancy, SMC: Seasonal malaria chemotherapy, IPTi: Intermittent preventive treatment in infants, PAAR: Priority above allocation request.
For each scenario, average trends in prevalence by microscopy, clinical incidence of uncomplicated malaria including both treated and untreated cases, and malaria-attributable mortality across all stochastic realizations were generated for U5 children and for all ages for years 2020-30. Case fatality rates were based on a review of published studies [75–78] and account for mortality due to severe malaria, maternal mortality from severe anemia in pregnancy, and mortality from malaria-attributable low birth weight in pregnancy (see Supplement). The 2020-30 projections for each scenario were compared against the simulated results for 2015 and against the 2020 BAU scenario to report projected relative difference in prevalence, uncomplicated malaria incidence, and malaria mortality among all ages and U5s for years 2025 and 2030.