Impact of intensified control on visceral leishmaniasis in a highly-endemic district of Bihar, India: an interrupted time series analysis

Visceral leishmaniasis (VL) is declining in India and the World Health Organization’s (WHO) 2020 ‘elimination as a public health problem’ target has nearly been achieved. Intensified combined interventions might help reach elimination, but their impact has not been assessed. WHO’s Neglected Tropical Diseases 2021–2030 roadmap provides an opportunity to revisit VL control strategies. We estimated the combined effect of a district-wide pilot of intensified interventions in the highly-endemic Vaishali district, where cases fell from 3,598 in 2012–2014 to 762 in 2015–2017. The intensified control approach comprised indoor residual spraying with improved supervision; VL-specific training for accredited social health activists to reduce onset-to-diagnosis time; and increased Information Education & Communication activities in the community. We compared the rate of incidence decrease in Vaishali to other districts in Bihar state via an interrupted time series analysis with a spatiotemporal model informed by previous VL epidemiological estimates. Changes in Vaishali’s rank among Bihar’s endemic districts in terms of monthly incidence showed a change pre-pilot (3rd highest out of 33 reporting districts) vs. during the pilot (9th) (p<1e-10). The rate of decline in Vaishali’s incidence saw no change in rank at 11th highest, both pre-pilot & during the pilot. Counterfactual model simulations suggest an estimated median of 352 cases (IQR 234–477) were averted by the Vaishali pilot between January 2015 and December 2017, which was robust to modest changes in the onset-to-diagnosis distribution. Strengthening control strategies may have precipitated a substantial change in VL incidence in Vaishali and suggests this approach should be piloted in other highly-endemic districts.

passed and only ∼2% of blocks are still above the target (February 2021), but resurgence may still occur-as seen in previous decadal cycles (Dye and Wolpert, 1988). Therefore, analysis from this recent intensified pilot is still relevant for future control policy.

Routine VL control in Bihar state
Current interventions implemented by the National Vector Borne Disease Control Programme (NVBDCP) involve biannual indoor residual spraying (IRS) of insecticide at state-level, passive case detection (PCD) at (block-level) by primary health centres and active case detection (ACD) by accredited social health activists (ASHAs); or via annual mobile camps (Supplementary information (SI) §S1).

Existing research base for interventions
An IRS review in Bangladesh, Nepal & India showed IRS had an impact on sandfly densities when properly conducted but did not show a significant impact on VL case incidence (Picado et al., 2012). The only large-scale randomised control trial of a vector control intervention on infection incidence (the KALANET project) found no evidence that large-scale distribution of long-lasting insecticidal nets provided additional protection over existing control practices (Picado et al., 2010(Picado et al., , 2015. A multi-site ACD screening intervention by ASHAs in highlyendemic Muzaffarpur & Saran districts discovered 6·7-17·1% more cases than PCD alone (Hirve et al., 2010). Overall, robust evidence is lacking on intervention effectiveness from field trials. Nevertheless, we hypothesise that a combination of strengthening the ACD referral system through VL-specific training for ASHAs, higher quality IRS by well-trained & supervised spray teams, and information, education & communication (IEC) community activities (SI §S2) could produce measurable incidence reductions. It is expensive to run a control programme of this scale: requiring coordination between Rajendra Memorial Research Institute of Medical Sciences (RMRIMS) and the Ministry of Health & Family Welfare, for administrative & logistical support and needing 166 spray squads for several weeks, twice a year (SI §S2). Therefore, any policy decision to apply this costly intervention to highly-endemic districts in future requires an evidence base.

Intensified control in Vaishali district
RMRIMS conducted an observational study on the impact of intensified VL control covering 1,569 villages in all 16 blocks in Vaishali district in late 2014-early 2015 (when 15 blocks were above the elimination threshold), while standard control by the NVBDCP continued in other districts ( Fig. 1 & SI §S1-2) (Kumar et al., 2020). The triad of ongoing interventions, which began asynchronously, are specialised ASHA training

Research questions & rationale for spatiotemporal model
In this study we estimate: (i) whether intensified control additionally contributed to the decline in VL cases in Vaishali, versus other districts (RQ1), & (ii) how many VL cases were averted by the pilot? (RQ2). Answering these questions is complicated since incidence was already falling in Vaishali before the pilot started (Fig. 2). Crude calculations indicate decreasing case counts year-on-year: 664 in 2014, falling by 38·1% to 411 in 2015, and by 56·4% to 179 in 2016 (Kumar et al., 2017). To estimate the impact of the pilot while accounting for the decreasing background secular trend, we compared Vaishali with other districts rather than analysing it in isolation. The model is informed by prior VL epidemiology & spatiotemporal features of the setting (SI §S6-7) (Bern et al., 2010). To estimate the number of cases averted, we fit the same model to a subset of pre-intervention months and make counterfactual predictions of case counts with which observed case counts can be compared (Graphical abstract & SI §S11).

Longitudinal dataset
Monthly VL case counts (by diagnosis date) for 33 out of the 38 districts of Bihar from January 2012-December 2017 were provided by the State Vector Borne Disease Office (Government of India, 2017). Our analysis included HIV-VL cases from January 2015-December 2016 and HIV/TB-VL cases from January 2017-December 2017 but excluded post-kala-azar dermal leishmaniasis (PKDL) cases (SI Fig. S9). The 33 study districts formed a contiguous island of transmission without the five remaining districts (Aurangabad, Gaya, Jamui, Kaimur & Rohtas, which are considered non-endemic) (Fig. 1a). Monthly district populations were estimated from 2001& 2011censuses (Government of India, 2015. District shapefiles provided adjacency information (GADM, 2015). The Institutional Ethical Committee of RMRIMS approved the intensified control programme (03/RMRI/EC/2018). University of Warwick's Biomedical & Scientific Research Ethics Committee (REGO-2018-2231) approved this analysis.

Descriptive analysis
Districts were compared by their ranked incidence levels and yearon-year changes in monthly incidence (SI §S3). Changes in rank position enabled us to crudely compare the relative changes of Vaishali to other districts, in the context of a state-wide medium-term decline in incidence. Using the two-sample two-tailed Wilcoxon test with continuity correction, we assessed if the ranks before & during the pilot were different. Evidence for global spatial correlation in incidence was assessed before & during the pilot with a Global Moran's I statistic hypothesis test (SI §S4). The effective reproduction number̂( ) for Vaishali & non-pilot districts was estimated to explore temporal patterns in transmission that may have been affected by interventions or seasonality (SI §S5) (Cori et al., 2013(Cori et al., , 2020.

Interrupted time series analysis (ITSA)
ITSA is a subset of regression discontinuity analysis which we applied to districts' longitudinal case counts to assess the impact of this non-randomised pilot while adjusting for existing trends (Kontopantelis et al., 2015;Bärnighausen et al., 2017;Bernal et al., 2017). The assignment variable in this ITSA is the calendar time (Eqn. 2:2) of the start of the pilot implementation in Vaishali. We evaluated the dynamics of case counts before & during the pilot using a spatiotemporal framework (Held et al., 2005Meyer et al., 2017;Bracher and Held, 2020;.
The process producing observed cases , in any district in diagnosis month is assumed to follow a Negative Binomial distribution with mean , & variance 2 , conditional on a weighted sum of cases from the previous 12 months ∑ 12

=1
, − , where is the weight for the cases months ago (i.e. distributed-lag autoregression) (Eqn. 1 & SI §S7). This distributed-lag distribution represents the diagnosis-todiagnosis (DD) distribution, i.e. the distribution of times between VL diagnoses of infector & infectee (SI §S6), akin to a 'diagnosis' serial interval distribution. It better represents the temporal correlation of diagnosis times than a naïve lag-1 autoregression (Bracher and Held, 2020). The normalised DD distribution is informed by an estimated incubation period (mean = 6 mo) (Chapman et al., 2015, which broadly agrees with literature estimates (Boelaert and Sundar, 2014), and an onset-to-diagnosis (OD) distribution (mean = 1·47 mo) from a Bihar study in the third quarter of 2012 (Jervis et al., 2017). With the DD distribution being a central assumption of all our models, we focussed a sensitivity analysis on the OD distribution to assess impact on RQ1 & RQ2 (SI §S12).
We expanded the base model with the pilot effect ( apriori primary variable, ( ) pilot ) and annual seasonality in the epidemic & endemic components (Eqns. 2:1 & 2:2, SI §S9) to form the final model: , and corrections for the first 12 months of the pilot due to delayed-lag intervention effects (SI §S9).
Model selection was based on the lowest Akaike information criterion (AIC) of each candidate model versus the best-performing model from the previous selection step, while monitoring changes in parameter uncertainty, particularly for the primary variable (Paul and Held, To assess the sensitivity of the final model's parameters to , we also reported their range when the start month was varied. For model validation (SI §S10) predictive performance in the prepilot period was visually assessed through fanplots & probability integral transform (PIT) histograms (Czado et al., 2009;Paul and Held, 2011;Abel, 2015;Meyer et al., 2017). The fitted values of the final model's case counts were plotted against their residuals to assess heteroskedasticity.

Counterfactual model
A counterfactual model, formed by omitting from Eqn. 2:2 the ( term, was used to predict the number of cases that would have occurred had there been no intensified control in Vaishali. Informed by the most likely start month * , we estimated the cases averted in Vaishali during January 2015-December 2017 by summing the monthly differences between simulated case counts from this counterfactual (fitted on January 2013-December 2014), and the final model (SI §S11). Cases averted were also presented as a percentage of those that would have occurred under the counterfactual model (Graphical abstract:D).
When developing a single model to both infer the pilot effect (RQ1) & generate counterfactual predictions (RQ2), it is unclear how to weight AIC & predictive performance (SI §S10) for these respective purposes (Nightingale et al., 2020). So pragmatically, we optimised for fit first; and prediction second, working from the pilot model. For model validation, the base, pilot & counterfactual models' goodness of fit were compared by AIC for the period January 2013-December 2017 (SI §S10). The predictive performance of the base & final models for January 2015-November 2017 was also compared (SI §S9).

Trends in diagnoses
Case counts fell in most districts, including Vaishali, from 2012-2017 (Fig. 2). During the pilot years 2015-2017, monthly cases in Vaishali declined substantially in absolute terms compared with cases in its highly-endemic neighbours and the district mean of the rest of the state (Figs. 1a & 2). Both during & before the pilot, Vaishali had the 11th largest yearon-year percentage reduction in monthly VL incidence out of 33 reporting districts (averaged over 36 monthly incidence ranks from 2015-2017 and over 24 months from 2013-2014 respectively). It also had the 9th highest VL incidence during the pilot period (averaged over 36 monthly incidence ranks from 2015-2017) versus pre-pilot when it was the 3 rd highest (averaged over 36 months from 2012-2014, < 1e-10 Wilcoxon test).

Seasonality & spatial correlation
Across Bihar, an annual seasonality in case counts was apparent, whose signal weakened as endemicity fell (Fig. 2). However, at district-level the strength of the seasonal signal varied and was only recognisable for some high-endemicity districts (e.g. Saran had a strong seasonal signal while others did not-note not presented in figures.).
Spatial correlation in incidence between neighbouring districts was apparent both before & during the pilot; Global Moran's = 0·36, = 0·002 (10,000 simulations) and = 0·40, = 8e-04, respectively. This supports the use of the between-district neighbourhood component in Eqn. 2:3. Vaishali was surrounded by neighbours with a range of endemicities, which either remained constant (e.g. Saran) or declined. Although incidence was declining in Vaishali, it was also declining in many other districts, yet clustering remained among the other districts, while Vaishali was dissimilar to its neighbours.

Effective reproduction number
The estimated district-specific effective reproduction numberŝ( ) generally follow an annual seasonality (Fig. 3b) which supports using seasonality in the model (Eqns. 2:1 & 2:2). Compared to the average trend of the other 32 districts, Vaishali saw sustained̂< 1 during 2012/3, with the second noticeable sustained reduction around the pilot start (Fig. 3a): after summer 2015, Vaishali'ŝdid not return to a seasonal peak around January 2016 unlike the mean of the other 32 districts. However, this effect only lasted the 2015/6 season and Vaishali'ŝresurged at the end of 2016. Given 2017's lower incidence, the impact of this above-onêin terms of new cases would have been less than if it had occurred at 2015's case levels.

Pilot model estimation
The final model selected consisted of a Negative Binomial distribution with population offset, annual sinusoid in epidemic & endemic components to account for seasonality, time-specific endemic intercept for high/low incidence, a distributed-lag epidemic component with a single change-of-intercept in Vaishali in January 2015, a constant distributed-lag contribution from directly-adjacent districts in the neighbourhood component, fixed intercept means in the epidemic component (one for Vaishali & one for the other 32 districts), and overdispersion by high/low-endemicity districts (Eqns. 1 & 2, SI §S6-9 & Fig. S4). The final model fitted better than the base model ( AIC = −308·3) and showed a reasonable fit to the observed case counts for Vaishali (SI Fig. S5) but across all districts was prone to overestimating low counts (SI Fig. S6). We chose January 2015 as the pilot start month * as it had the lowest AIC. Table 1 shows the parameter estimates for the intervention effect which can be interpreted as follows. For Vaishali pre-pilot, an estimated average 67·8% of the weighted sum of the previous 12 months' case counts contributed towards the current month's case count, versus 70·8% (95%CI 67·7-73·8%) for other districts. This means a hypothetical same-sized epidemic in any of the other districts would take slightly longer to die out on average than if it was to occur in Vaishali. During the pilot, this 67·8% contribution was estimated to fall (by 27·3%) to 49·3% for January 2015 onwards, where the pilot effect confidence interval represents a significant drop (i.e. 8·8-45·8%). The . The inferred infection times on the x axes was calculated by subtracting the mean incubation period and mean onset-to-diagnosis time (7-month total shift) from the diagnosis month. The dashed vertical line indicates the length of the 7-month sliding window used to smooth thêestimates: before this date estimates are unreliable as they only include partial data within this sliding interval. The dotted vertical line indicates the start of the modelled intervention. estimated endemic contribution per district since January 2012 (based on a mean district population of 2·8 million) was practically nil (∼0 cases/mo) for low-incidence settings (<11 observed cases/mo), whereas high-incidence settings (≥11 observed cases/mo) are estimated to get 3/6 cases/mo for low/high seasons, respectively. In absolute terms the seasonality term contributed more to the epidemic component than endemic (SI Fig. S5)-0·77 times at the November minimum & 1·31 times at the May maximum of the epidemic component. Each district received an estimated 0·5% average contribution from each of the adjacent districts' weighted sum of their previous 12 months' cases. The standard errors of all parameters were within reasonable bounds. Parameters were mostly insensitive to pilot start month , however, the pilot effect on the epidemic component ( ) pilot and Vaishalispecific intercept ( ) Vaishali did differ by up to 13% & 7%, respectively, versus their value for a January 2015 start. The changes in the mean parameter of OD distribution had negligible effect on both RQ1 & RQ2 (SI Table S2).

Estimating cases averted
The counterfactual model showed reasonable predictive performance (SI §S9-10) before the pilot, notwithstanding the last four months of 2014, where the limited duration of this test period caused convergence issues (SI Fig. S7). The final model (SI Fig. S8b) produced forward predictions for 2016-2017 generally sharper & closer to the observed time series than those of the counterfactual model (SI Fig. S8a). Predictions of both models into 2017 were less robust to Kumar, Siddiqui, Pollington et al.  large departures from the mean trend, e.g. the unexpected July 2017 peak in observed cases (SI Fig. S8b), as the epidemic component was diminished by this point. Predictive performance of the final model was poorer at extrema, especially high counts (SI Fig. S10).
Simulations comparing pilot & counterfactual models suggest a median 352 (IQR 234-477) cases were averted in Vaishali during the pilot from January 2015 (2% of 100,000 simulations had estimated negative cases averted), which would have accounted for an estimated 31% of cases if there had been no intensified control (Fig. 4).

Discussion
This study comes at a critical point in VL elimination, where highendemicity districts are predicted to be the hardest in which to reach the elimination target (Le Rutte et al., 2018). Our analysis of the Vaishali pilot study suggests that combining existing interventions with special attention to quality, might contribute to additional reductions in VL incidence.
Descriptive analysis suggests a significant change in the case counts in Vaishali for the first two pilot years 2015-2016 relative to other districts, which is supported by our detailed spatiotemporal analysis that accounts for decreasing trends in cases pre-pilot and neighbouring district effects. When the study started, 15 out of 16 blocks in Vaishali were above the elimination target of 1 case/10,000 people/year, but all blocks apart from Raghopur (where flooding interrupted the pilot in August 2017) are now below the target. Model simulations characterising the pilot period suggest that several hundred cases have been averted since 2015, which was robust to changes in the OD distribution.
We cannot conclusively attribute the additional decline in case counts in Vaishali from 2015 to the intensified control programme because this is an observational study. For internal validity of an ITSA, the continuity assumption must be met so that one is reasonably confident that ''no other interventions or confounding covariates than the treatment of interest in analyses changed'' at the intervention start month (Bärnighausen et al., 2017). As the pilot & initial decline were concurrent and because no other widespread interventions were in place (SI §S1), we conclude that the additional decline was most likely due to the intensified interventions.

Limitations
This study does not apportion how much each of the pilot's triad of interventions contributed to the decline nor does it include covariates that describe the time-varying susceptibility of sandflies to the deployed insecticides. Modelling suggests this pilot's high 90% household coverage per block (SI §S2:2) would have been insufficient alone to reach disease elimination (Fortunato et al., 2021). In addition, a recent study in two highly-endemic districts of Bihar suggests IRS, as implemented under the national control programme, has a negligible impact on sandfly abundance (Poché et al., 2018). Vector abundance, insecticide susceptibility & IRS coverage data from Bihar's districts would allow further investigation.
'Single-world' matching of counterfactual simulations to their corresponding pilot simulations could produce a similar point estimate for cases averted but with lower stochastic variation (Kaminsky et al., 2019), producing the averted estimate with a narrower uncertainty band (c.f. Fig. 4) but is beyond this study's scope. A control group is also lacking as the 32 comparison districts could have unobserved confounders distributed heterogeneously across them, which limits the external validity of the analysis; in addition to inferences coming from a single district.
We do not know the treatment information of some Vaishali cases that chose nearby district hospitals nor other districts' cases migrating into Vaishali, which could affect the estimated contributions of the epidemic & neighbourhood terms in the model to Vaishali's case counts. It is also unclear how drug supply may have impacted incidence since the national programme introduced single-dose liposomal amphotericin B in 2015-2016. Some of the largest differences among the 32 non-Vaishali districts are the VL endemicity & mean OD (Jervis et al., 2017), however, our model does not account for these heterogeneities. If ASHA training reduced onset-to-diagnosis (OD) times, and thus infectious durations & subsequent incidence, this would have also shortened the DD distribution, meaning that our inferences & prediction are biased. However, we expect any large reductions in the infectious duration would only marginally affect the OD distribution as the mean infectious period was only 11% of the mean DD.
Underreporting of cases, estimated at 15%-18% in Vaishali in 2012-2013, with a non-uniform age distribution, may have affected our results (Das et al., 2016;Jervis et al., 2017). However, NVBDCP introduced mandatory VL reporting state-wide for the public sector on 7 January 2016. Although HIV/TB-VL coinfection data is included in the monthly cases, we have been unable to stratify their status in the model due to this data only being available since 2015. In a Vaishali district hospital in 2011-2013, VL admissions who were unknowingly HIV+, had OD times on average 3 weeks longer (Burza et al., 2014a); their underdiagnosed HIV-VL status accounted for 2·4% of admissions, rising to 5% in middle-aged men. This may also be important if HIV-VL-coinfected individuals contribute disproportionately to transmission (Burza et al., 2014b). If they do, then the pilot effect in 2017 for Vaishali, a district with a rising proportion of HIV-VL coinfections, may be underestimated. Furthermore, PKDL cases are not incorporated into the analysis as case counts were unavailable from 2012 but recent studies suggest they contribute significantly to transmission as VL incidence declines (Mondal et al., 2019;Chapman et al., 2020).
Despite these limitations, we recommend further pilots in highlyendemic settings with additional collection of time-varying district covariates and assessment of cost effectiveness. Commendably, the WHO roadmap now has three 2030 targets (WHO, 2020a,b): 1. PKDL elimination (all PKDL cases detected & treated from recovered VL cases followed up for 3 years) 2. VL case fatality rate (<1% nationally) 3. reaffirming the VL elimination target (all blocks in India at <1 new/relapsing VL case per 10,000 population).
Thus widening research questions to these recent goals would support India's efforts towards VL elimination.

Conclusion
Can intensified control reduce VL incidence more quickly in a highly-endemic district? Our robust analysis shows that observed VL case counts did fall more quickly in Vaishali district than other districts, in line with previous crude analyses (Kumar et al., 2017(Kumar et al., , 2020 and estimates an additional outcome indicator as 'cases averted'. Since the design of this study (2014), VL policy now covers PKDL burden & VL mortality. We believe there is justification for piloting our approach in other highly-endemic settings, contingent on improvements in study design & analysis ( § 4.1 & SI §S13), to meet these policy advances.

Declaration of competing interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:  Meyer et al., 2017), and two anonymous reviewers for valuable input. We also credit the open-source software we used (Pebesma and Bivand, 2005;Wickham, 2007;Plate and Heiberger, 2016;Wickham, 2016;Bivand and Wong, 2018;Meschiari, 2021;Garbuszus and Jeworutzki, 2021;RStudio, 2021;Bivand and Rundel, 2021 should not be attributed to these funders or any person connected with them. No funders had a role in the study design, collection, analysis, interpretation of data, writing of the report, or decision to submit the manuscript for publication. No others were paid to write this article. Final approval for publication: All authors. Kumar, Siddiqui, Pollington et al.

Appendix A. Supplementary data
Supplementary material related to this article can be found online at https://doi.org/10.1016/j.epidem.2022.100562. We detail the study design, analysis & modelling in SI §S1-14 (also found at DOI: 10.5281/ zenodo.5579140) and provide a policy-relevant items for reporting models in epidemiology of neglected tropical diseases (PRIME-NTD) summary (SI Table S3), to assist in communicating ''the quality & relevance of modelling to stakeholders'' (Behrend et al., 2020). Analysis code is available at github.com/t-pollington/ITSA.